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PREFACE 


This manual describes the mathematical routines of the FORTRAN Common Library which is part of 
FORTRAN Extended Version 4. It is assumed that the reader is familiar with FORTRAN Extended. 
FORTRAN Extended operates under control of the following operating systems: 

NOS 1 for the CONTROL DATA® CYBER 170 Models 171, 172, 173, 174, 175; CYBER 70 
Models 71, 72, 73, 74; and 6000 Series Computer Sytems. 

NOS/BE 1 for the CDC® CYBER 170 Series; CYBER 70 Models 71, 72, 73, 74; and 6000 
Series Computer Systems. 

SCOPE 2 for the CDC CYBER 170 Model 176; CYBER 70 Model 76; and 7600 Computer 
Systems. 


CDC manuals can be ordered from Control Data Literature and Distribution 
Services, 8100 East Bloomington Freeway, Minneapolis, Minnesota 55420. 


This product is intended for use only as described in this document. Control 
Data cannot be responsible for the proper functioning of undescribed features 
or parameters. 
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The Math Library concerns itself with computations upon four 
different number typest integer, single [precision floating-point], 
double [precision floating-point!, and complex [floating-point]. For 
each nymber type there is a well defined set of valid forms 
Treoresentations], each one representing a particular point on the real 
line or in the complex plane. In adcition, for each of the floating¬ 
point forms, there is a well-defined set of semi-valid forms, none of 
which represent numbers, but which instead give some indication of the 
nature of the [erroneous! computational process that produced them. ftll 
other bit configurations in words thought to contain numbers of some one 
of these types are termed invaI id . 

For these four number types, the valid, semi-valid, and invalid forms 

are 

1. Integer. 


Validl 


The ordinary, ore word, right-}ustified, one*s- 
complement binary represen tations of all integers 
from to ?*►«-! . Zero may be represented as 

either positive zero [all zero bits], or negative 
zero [all one bits!. 


Semi-valid* None 

Invalids Any bit configuration wherein the top 12 bits are not 

a fI the same. 


2. Single. 
Valid: 


The normalized, one word, forms 
Tloating-point representations, 
Computer Systems* Manual.) Zero ma 
as either positive zero, or negativ 


of the internal 
(See corresponding 
y be represented 
e zero. 
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Sem L-va I id 


The four forms known as positive infinite, negative 
infinite, positive indefinite, and negative 
inde finite . 

Invalid* ftny non-zero and non-semi-va!id bit configuration 
wherein bit and bit 59 are the same, etc. 


Double. 

Valid* The forms of the irternat, double-precision floating- 

Doint reoresentations wherein the first word is 
normalized and the second word either has an exponent 
that is 48 smaller than the first word, or, if that 
underflows is zero. The signs of both words must be 
the same except when the lower part underflows to 
zero. Zero may be represented as either positive 
zero or negative zero. 

5>emi-valid* The forms wherein the first word is a single semi- 
valid form. The second word may be anything. 

Invalid* Sign disagreement between the two words, first word 

an invalid single, second word with an exponent not 
as defined above, etc. 


4. Complex. 

Valid* All the two-word forms wherein each word is a valid 

single number. 

Semi-valid* All the two-word forms wherein one word is a semi- 

valid single number, and the other is either a valid 
or a semi-valid single number. 

Invalid* All the two-word forms wherein either word is an 

invalid single number. 

The following two general rules apply to the use of these number 
forms in computational operations, either within the Math Library or 
within FOI?TRAM compiled code 

1. Unless specially documented otherwise, if a valid form of the 
appropriate number type is employed in a computational operation, a 
valid number of the appropriate type will result. The documented 
exceptions to this cover such things as computing an answer which 
exceeds the limits of the valid forms, or performing a 
mathematically invalid operation. 

?, Unless specifically documented otherwise, if either* 

a. a semi-valid or invalid number is employed in a computational 
operation, or 

I 2 
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b. the documented limits in rule 1 above are exceeded^ 
then the result is undefined Ci»e. the program may continue without 
warnino, it may terminate abnormally with or without diagnostic* 
it may continue for a short period and then terminate* etc.l. 
The documented exceptions to this cover some cases wherein 
certain forms of checking are done, and also some cases wherein 
certain semi-valid forms are produced, etc. 

These two rules define the limits of COC support in the area, and 
also the completeness of COC supporting documentation. When a result is 
undefined, there is Qa guarantee that the actual behavior will be the 
same from run to run, or that it will remain constant under normal 
product maintenance. 


ITT. CLASSIFICATION OF ROUTINES AND CALLS 


The FORTRAN Common Library mathematical routines (abbreviated * math 
library routines) compute those mathematical functions explicitly 
mentioned in FOPT^AN. These functions may be divided into two classes - 
the intrinsic functions and the external functions. Intrinsic functions 
are simpler functions whose use speeds execution of programs and saves 
coding effort by occasioning replacement of freouently used sequences of 
F0*^TRAN statements with efficient in-line code during an intermediate 
assembly, or with calls by name to routines (when in traceback mode). 
The list of intrinsic functions appears in the Appendix. External 
functions are mathematically more sophisticated functions whose routines 
require more memory space and execution time. Calls to math library 
routines may be of two forms - calls by name and calls by value. When a 
routine is called by name, a parameter list is formed in memory, and the 
first-word-address of this list is entered into register A1 before a 
return fump is made to the routine. When the routine is called by 
value, the arguments are entered directly into operand registers 
according to certain rules, before a return Jump is made to 
the routine. The first word of the first argument is entered into XI, 

the first word of the second argument is entered into X3 and the first 

word of the third argument in X5. If an argument should be double¬ 
precision or complex and hence take two words, the second word is 
entered into the next register, viz. XE or X4. Lastly, the first word 
of a complex argument is always the real part, and the first word of a 
double-precision argument is the upper half. For calls by name and 
calls by value, the result of the computation is returned in registers 

X6 and X7, a one-word result being returned in X6, and the second-word 

of a two word result being returned in X71 • (A Juxtaposition symbol 
(A) is sometimes used in the documentation * it denotes that a two-word 
result occupies the two registers in the order indicated.) 


IV. TERMINOLOGY 


Some conventions have been introduced in this documentation. 
Symbolic names are always delimited by blanks, and any latin letters 
appearing therein are in upper case. a denotes Iuxtaposition . and is 


60498200 C 


3 


used in referring to coffolex or double-precision quantities. AH values 
given are in decimal, unless otherwise noted. Error shall meani 
(computed value - true value}. Rel ative error shall meant lerror/true 

value). An acau ro ent.set is an ordered n-tuple of arguments (xl,...,xn)y 

where xl,...,xn are the arguments in order. For convenience, we 
identify arguments with corresponding 1-member argument sets. The inout 
of a routine is the collection of all argument sets for which that 
routine has been designed to return a result meaningful to the user. 
For example, the inout range to SIN is the col Iection of all floating¬ 
point quantities, whereas the input range to SINCOS= at entry point SIN, 
is the collection of definite in-range floating-point quantities not 
exceeding pi.?**® in absolute value. PQS. INF. abbreviates 
7777,nqnO,0000,0000,OOOOB , NEG.INF. 4000,0000,0000,0000,0008 , 

abbreviates 1777, 0000,0000,0000,00008 , and NFG. TNOFF . 
abbreviates 8000,0000, 0000,0000,00003 . In this document, " routi n g '" 
shall mean the source code or the object code obtained from programs in 
the UPDATE library mentioned at the begining of this Introduction, 


V. EPROR GRAPHS 


The errors of a routine are composed of two partst the algorithm 
error, including errors in the coefficients used in the algorithm? and 
machine round-off errors. A curve representing the error due to the 
algorithm and its coefficients is usually a smooth, wavy curve with 
discontinuities at breaks in the range reduction technique. The error 
of the coefficients involved in range reduction may also show up. 
Usually, a good algorithm with good coefficients will not have an error 
bigger than one-half in the last bit of the result. Round-off (and/or 
trunca^^ion) error is difficult to predict or graph. Suppose f(x) were 
aoproximated by x+c^x® and x>>c®x®^ • Then by analyzing how an add 
instruction would work on x and c^x* , one finds that a few bits are 
dropped off after the last bit in the result. If rounded add is used 
then the resulting error is between -1/2 ard 1/2 in the last bit? the 
error in computing c’^x® makes it even worse, A graph of round-off error 
is so discontinuous that little can be done other than showing the 
maximum and minimum error over small intervals. 

The magritude of a relative error can be analyzed in two ways* 
relative error = (routine - ex act)/exact? or figuring out how many bits 
the routine differs from the exact value ("bit error”). In the first 
case, we are talking about single precision algorithms accurate to less 
than 2F-15 (usually) and rourd-off errors less than lOE-15 (usually). 
Note* changing the last bit in a single precision number produces a 
relative change of between 3,5E-15 (for a large mantissa) and 7.16-15 
(for a small, but still normalized, mantissa). In determining how many 
bits off a routine is, the function is evaluated in double precision and 
this is rounded to single? then (assuming the exponents are the same) 
the mantissas are subtracted and the integer difference is the bit 
error. 
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a typical plot covers ore s ingl e*argumer.t» single-precision function 
over a range of argument values (plotted linearly or logarithmically) 
with the ordinate ranging from -llE-15 to llE-15 representing relative 
error. The saw-tooth curves represent places at which relative error is 
f t !V2 « and 1/2 bit error. Discontinuities occur where the 

routine produces a result that is a power of 2 ; the argument values are 
given (they are found empirically* so only an appropriate number of 
digits is printed). 

any point that is between the -1/2 and 1/2 saw-tooth curves 
represents a case of the routine being as accurate as possible* anything 
between 1/2 and ?/2 is 1 bit high? etc. 

(In algorithm error curve wiggles around through the middle of the 
plot, Tt shows the relative error of the algorithm over the given 
argument range. Its discontinuities are usually due to the range 
reduction part of the algorithm. For this curve* the algorithm error is 
(alg - exacty/exact where alg is routine rewritten to use double 
precision operators instead of single but keening single precision 
coefficients single. Therefore it incorporates such things ass a 
polynomial can't auite equal a transcendental function and oi/2 can't be 
represented exactly. The coordinates of the highest point are indicated 
next to it. 

The overall error is bounded (empirically) hy two lagged curves with 
arrowheads on them. The number of different arguments fed to the 
function is given on the olot? each corresponding point is either at the 
tio of one of the arrowheads or strictly between the pair of curves. It 
is quite oossible* even likely* that there are points which do not lie 
between the two curves. However* one could* with reservations* assume 
♦he curves are "close** to true least upper bound and greatest lower 
bound curves. 

The arguments are chosen randomly as follows, ftfter starting with 
the smallest argument, each argument is the previous argument olus 
P/\NF{n)»k * where k is a constant. On a logarithmic scale this 
algorithm is approoriateIy modified so as to get an even distribution on 
the resulting olot. 

Note that "ordinary** numbers (rational numbers, multiples of log 2 or 
pi, etc.) probably will not be sampled. 

There are usually about 25(5 points (arrowheads) on each of the 
houndlno curves. The algorithm for finding arrowheads goes as follows. 
Given arrowheads x and y , the last two on the list, point z (formed by 
an argument and the relative error of the routine for that value) is 
added to the arrowhead list if xyz forms a convex curve or the abscissa 
of X and z are "too far" apart. Otherwise* arrowhead y is deleted from 
♦he list and the test for inclusion is retried. Points going beyond 
11^-15 are forced to the boundary. The largest relative error 
encountered is labeled with its coordinates. Various statistics are 
printed concerning the distribution of points. The percent within each 
bin of width lE-15 with the percent above lOE-15 (below -lOE-15) being 
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stated between lOE-15 and llE-15 <-llE-15 and -lOE-151. Sit errors are 
similarly handled (with anything above 3 being put with 3). Eiapty bins 
are not listed The "HEAN R«E»“ is the mean of all ordinates* The "RMS 
P.F." is 

STRTCIsum of PE**2)/(number of pointsi - (MEAN R*E.)**2 « 

i.e., the standard deviation of relative error* 


Hoji ,LQ-.2aad-a-.£lal 

Here are some cause-and-effect statements* by taking the inverse of 
the statement one has a way to look at a plot and deduce what the 
algorithm is doing. 

1. If f(x) = 2**n * (x*g(x» where g(xl is small compared to x and 
rounded add is used« then the bounding curves will roughly 
paralled the algorithm error and will be as far apart as the 
inner saw-tooth curves* (Unrounded add would transpose the 
curves by 1/2 bit.) 

2. If f(x) = ceg(x) then the bounds will be transposed by the error 
in c. 

3. If f(x) = c»g(x) then the distance between the bounds for f(x) 
will usually be wider than for g(x) I in particular f(x) will 
probably have bounds at least 2 bits apart. 

4. If f(x) = g(x) <-(h(x) + d(x) ) where g, h, or d may be constant and 
one of the additions produces ar unnormalized result, then the 
bound curves may be translated and/or spread farther apart than 
for a nearby area where the addition happens to be normalized. 

E. If f(x) is broken into numerous sub-irtervals (e.g* 16), then 
the algorithm error curve will be dominated by discontinuous 
Jumps in the constants used for table lookup* 

VI. MISCELLANEOUS FACTS 


Arguments of trigonometric functions and results of inverse 
trigonometric functions are always measured in radians. Some statistics 
concerning the UPDATE library of mathematical routines are given for the 
Cytfr 74 , There are 128 routines. The central memory required to 
UPDATE all routines is 3630!) (octal) words. The central memory required 
to assemble all routines is 50700 (octal) words. The time required on a 
r;Y;^c-g 74 to UPDATE all routines is 4.8 CP seconds, and the time required 
to assemble all routines under COMPASS is 28 CP seconds. The average 
assembly time for individual routines is .2? CP seconds. These times 
wilI be shorter on the CYBER 76 and longer on the CYBER 72 and 73. 
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’ROUTINE'S EUNCTTON. 


1 . 1 . 


Type. FORTRAN external 
floating-point argument. 


functions. The routine accepts a 
and returns a floating-point result. 


1 . 2 . 


Purpose. To accept calls from FTN compiled code for 
computation of the inverse cosine and inverse sine functions. 


2. MFTHOn. 


The input range is the collection of all valid floating-point 
quantities in the interval f-1.,1.1. Arguments outside this range 
will initiate error processing. 

Formulae used in the routine arel 


arcs in(x) 
arcos(x) 
arcsin(II 
arcos111 
arcsin 1x1 
arcos (xl 


= -arcsinC-x) 

= oi-arcos(-xl 
- pi/2 
= 0 

= pi/2-arcos(xl 
= arcos(l-g(x,nl|/2**n 


x<-.5 

x<-.5 


.5<x<.l 
• 5<x51. 


where 


g(x,OI a 1-x 

q(x*n^ll a <fg(x«n) - 2g<x»nl*. 

arcos(x) a pi/2 - arcsinix) -.55x?.5 

arcsin(x) a x«'X^*s*( (wez-f l^weaetr/(e-x^l I 

-.5<x<.F 


where 


w js Ix^-cl^z+'k 


and 

z a (x^erlX^ei 

The constants employed are* 

r a T.17317007853713 
e = 1.16039462573gn2 
m = 59.3190559607983 
c a -2.36958885561288 
i = 8.22646797079917 
1 a -39,6294815974555 

4 = 37.4592309257582 
a = 349.319357025144 

5 = .746926199735419 * 10**-3 
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The approximation to arcsln is an economized approximation 

within the class obtained by varying r«e«mt*«.»s« The algorithm 
employed is as follows* The argument x is supplied to ACOS* or 
ASIN* in XI « and the result is returned in X6 • 

a* If ftCOS* entry* go to step g* 

b* If lxl>*5* go to step h* 

c. n«‘0 fLoop counterl 

u^-O if AS IN* entry 
••pi/2 if ACOS* entry 

d* z*-lverl ’yti 
w«*Cvc> •z«-K 

p«-q*s*q*y*f fw*z-J I mweatiiyCe-y n 

p*>u-p 

Xfi*-p/2**n 

e* If ASIN* entry* go to step k* 

f* If X is in f'-*5*l*l* return* 

X6«-2*u-(X6» 

Return* 

g* If Ixl < *5* go to step c* 

h* If X = ♦l*-l or X is invalid* go to step t* 
n'-O (Loop counter! 
y’-l-lxl* and normalize y * 

i* h«-4*v-2*y* 
n*-n*l 

If 2*y<2-sqrt (31 * *267949192431* y^-h and go to step i* 

I* q«*l*h* and normalize q* 

y*-g» 

u*-pi/2 

Go to step d* 

H* XB^-u-IXfl* and normalize X6 * 

Affix sign of x to X6 • 

Return* 

I* If X #♦!* or -1* * go to step m* 

X6«*pi/2 If X = I* 

•-Pi/2 if X s -1* 

If ASIN* entry* return* 

X6«-0 if X = 1* 

►piifx*-!* 

Return* 
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m. Plug ACOS. entry point with ftSIN* entry point. If ASIN* entry. 
Initiate error processing. 

Return through ACOS. entry point. 


. ERROR ANALYSIS. 


The maximum absolute value of relative error of the approximation 
above to arcsin over I-*.5,.S1 Is 1.996*10**-15 . A graph of the 
relative error of this approximation is given in figure 6. Upper 
bounds on the absolute value of relative error due to machine error 
have been established in the following casest 


arcsin 

on 

(-.5,.5) 

- 

9.232 

* 

arc os 

on 

(- 5,.5) 

- 

1.67 3 

* 

arcsin 

on 

(-l.,l.l 

- 

4.050 

* 

arcos 

on 

(-1.,1.1 


1.618 

» 

The corresponding upper 
error in the routine are? 

bounds 

arcsin 

on 

(-.5,.5) 


1.123 

* 

arcos 

on 

(-.5,.5) 

- 

1.873 

* 

arcsin 

on 

(-1.,1.) 

- 

4.250 

* 

arcos 

on 

(-l.,l.l 

- 

1.638 

♦ 


IQ**-!** 

10**-13 / 

on the absolute value of relative 


10**.14 

10**-14 

1Q**.14 


For groups of 1000 arguments chosen randomly from the following 
intervals, the following statistics on relative error were observed. 


Entry 

Point 

Interval's 
Lower 

Bound 

In terval"s 
Upper 

Bound 

Mean 

Standard 

Deviation 

Minumum 

Maximum 

ACOS. 

-.5 

.5 

-9.435E-16 

1.547E-15 

-5.781E-15 

- 3.856E-15 


-1. 

-.5 

-4.331E-16 

1.746E-15 

-4.520E-15 

4.546E-15 


.5 

1. 

-5.098E-16 

1.843E-15 

-7.150E-15 

9.559E-15 

ASIN. 

" .5 

.5 

8.401E-16 

1.666E-15 

-5.328E-15 

4.916E-15 


-1. 

-.5 

6.209E-16 

3.268E-15 

-7.061E-15 

1.489E-14 


.5 

1. 

7.311E-16 

3.307E-15 

-7.160E-15 

1.554E-14 


6.1. ALGORITHM ERROR. 

For ASIN (x) , x in (-.5,.51 the error curve is depicted In 
the ASIN plot between 0 and .9 . (All of the ASIN plot is 
symetric about 0.) • The reason for it not being balanced 
around the axis is because the Chebyshev coefficient for x was 
thrown away and 1.0 implicitly used instead. For ASIN outside 
(-.5,.51 and for ACOS , t'here is range reduction first? this 
produces no algorithm error. At the end of the computation, 
some multiple of pi gets involved? hence, the curves are 
offset by an amount dependant on the error in pi. There are 
breaks in the algorithm error curve at plus/minus .5 , 
S0RT(3)/2 = .866025 , .9665926 , .991445 , .997859 , etc. 
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Each of these Is SQRT flU-orevious)/2J . 

6»2. Total Error. 

In flSIN for x in (-.5*.51 the routine boils down to x*-x3*(,,,> 
y hence the total error is dominated by that final addition 
and the total error curve closely follows the algorithm error 
curve plus/minus 1/2 bit. For x in C.5».866) the algorithm is 
as follows* y=l-x « z=<l-4v) «-2y* « ASIN<x)=pi/2*(oi/2- 
(z+7=*(...)))/2 . y is in C.5y.l34) , z is in (-.5,.5) . 
Nothing is lost in computing y y little is lost with z y and 
some is lost in the final part. The big lump when x is in 
(.5y.540Tn2l is caused by pi/2-(z+...) being greater than 
2. ? elsewhere it is less. This peak shows up at other places 
Cin ASIN noticeably in (.866y.878l and ACOS Just below each 
peak in the bit error curve) because of folding into (.5y.54). 
ASIN gets better near 1.0 because pi/2 predominates the final 
val ue. 

ACOS y except near 1.0 « is predominated by pi/2 . In 
particular* for x in f-l.y.5)y oi/2 is added on twice* first 
rounded then unrounded in order to give a near-perfect 
distribution. Near x=1.0 * so much folding goes on that a 
rather bad error is built up even before evaluating the 
polynomial. The graph gives an indication of the infrequency 
of error but does not show a worst case (15E-15 relative error 
has been experienced). 


4. EFFECT OF ARGUMENT ERROR, 


If a small error e occurs in the argument x, the error in the result 
is oiven approximately by e/Cl-x*)**.5 for ASIN and by -e/(1-x*)**.5 
for ACCS . 
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ROUTINE'S FUNCTION. 


1.1 


I,?., 


Type. A FORTRAN external 
point argument and returns a 


function. It accepts a floating- 
floating-point result. 


Purpose. To accept calls by name and by value for ALOG fat 
entry points ALOG . ALOGIO « ALOG. and ALOGIO, ) from FORTRAN 
programs. ALOG computes the natural logarithm function at 
entry points ALOG and ALOG. . and the common logarithm 
function at entry points ALOGIO and ALOGIO. . 


?. METHOD. 


The input 
in-range 


range to this 
non-negat iwe 


of all definite 
quantities. Upon 

i s 


routine is the collection 
non-zero floating-point 
entry, the argument x is out in the form x = y ♦ 2 n, where n 
an integer, and l.<y<2. . Then log x is evaluated from 
log X = log y ♦ 3/4*n ♦ (log 2 - 3'4)*n, 
where log y is evaluated as follows. The interval 11., 2.) 
divided up into the subintervals 
(1., i.l07?3fl769FTll♦ 

(1.107238769521, 1.357238769531), 

(1.607238769531, 1.857238769531), 

(1.357238769531, 1.607238769531), and (1.857238769531, 2.). 


IS 


'Centre 


poin ts' 


1.735100002271352, 

2. 


are chosen wi 

in suhinterval (a, 

b) 

with centre poin 

log y = log c 

where 

♦ 1 

oq 

((l*t)/(l-t)» 

t - (y - c)/( 

y ♦ 

c). 


log ((1 + t)/(I - 

t)i 

is 

then computed 

log «(l*t)/(l-t») 

. 

cv 

II 

♦t 

♦ C(3)*t3 ♦ C 


1.225803196513098, 

in these 
c, log y 


1.475803239208091, 
intervals. It y is 
is computed from 


by 

C(5)*t5 


♦ c(7)»tT 4 . c(9)»t» 


The 

the 


coefficients c(.3), c(5), c(7) and c(9) are chosen 
Taylor series for log ((l*t)7(1-t)) after the 


by truncating 
11th term, and 


taking a Chebyshev economization 


1arqest 

i 

nt er 

val symmetric about 

The con 

St 

an ts 

are 

c ( 

3) 

— • 

666666666666105 

c ( 

5) 

^ • 

4000000018947 

c ( 

7) 

— • 

2857120487 

c(9J 


22330022 


to a 9th degree polynomial over 
the origin which is applicable. 


the 


If the argument x is 
5ySAin= , and POS.INDEF. 


invalid, an error message is issued through 
is returned. 


ER»CP ANALYSIS. 


(We carry out the error analysis for computation of ALOG only. 
Bounds on machine error are the same for ALOG and ALOGIO here, while 
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the the graph of algorithm error for flLOGiO may be obtained from the 
graph for ALOG by multiplying by logfeJlO.) The maximum absolute 
value of the relative error in the algorithm over the interval II.» 
?.» is 1.698 * 10 ** -16. for entry points ALOG and ALOG. . The 
maximum absolute value error in the algorithm over the interval II.. 
?•> is 1.667 * 10 ** -17. A graph of the error in the algorithm 
over II.. ?,} is given in figure 8. An upper bound has been 
established for the absolute value of the error in the routine due 
to machine error at 5.045 ♦ 10 -14 * u. where u is the greatest 
integral power of 2. not exceeding the result. Hence an upper 
bound on the absolute value of the relative error in the routine is 
5.06? ♦ 10 -14. 

For groups of 10000 arguments chosen randomly from the following 
intervals at the entry points listed, the following statistics or 
relative error were observed. 


Entry 

Point 

In terval 
from to 

Mean’ 

Standard 

Deviation 

Mini mum 

Maximum 

ALOG. 

1 . 

.5 

.5 

.0001 

10*»-290 

2. 

2. 

1. 

1000. 

10322 

1.743E-16 

2.325E-16 

4.101E-17 

4.522E-16 

1.228E-15 

2.2a6E-15 

2.279E-15 

2.488E-i5 

2.223E-15 

1.439E-15 

-9.040E-15 

-1.058E-14 

-9.450E-15 

-5.562E-15 

-1.616E-15 

6.ig4E-15 
8.665E-15 
8.637E-15 
5.234E-15 
4.001E-15 

ALOG10 

.1. 

.5 

.5 

.0001 

10»*-290 

2. 

2. 

1 . 

1003. 

10322 

-2.726E-15 
-2.689E-15 
-2.826E-15 
-1.795E-15 
-2,Q15E-15 

2.723E-15 

2.770E-15 

2.897E-15 

2.526E-15 

2.1/8E-15 

-1.447E-14 

-1.346E-14 

-1.546E-14 

-9.208E-15 

-7.389E-15 

4.640E-15 

6.506E-15 

9.353E-15 

5.058E-15 

3.453E-15 

6.1. 

ALGORITHM 

ERROR. 






P ange 

re duct 

ion first 

to 1 ds 

arguments 

into 


(.9286194,1.85723911 the unfolding involves an approximate 
constant involving log 2 I hence. the error graph shows 
discrete lumps at 2*»r*l.857239 in the algorithm error plot. 
Further range reduction into the subintervals described above 
involves the use of log c. The values of c were chosen so 
that the 48-bit representation of log c would be correct to at 
least 59 bits. Hence, no noticeable error is caused by 
reducing into the subinterva 1s. Within each subinterval a 
polynomial is used* the polynomial is accurate enough to show 
essentially no error except near 1.107239 . 

6.2. TOTAL EPROR. 

The final computation is log x = (C<(att)+t)+pl+b)+b where 
a = log 2 - 3/4)*n . 
p = c(3>»t3*-.,. , and 
b = (3/4 » n log cl/? . 

In general p<t<a<b except that a and/or b could be zero. The 
order was chosen in order to minimize error accutrul at ion. b 
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is added in twice in order to cut down on error and eliminate 
a normalize. Because of all this adding going on» the error 
graph lumps around at odd times and by fairly small amounts, 
ca lump probably corresponds to a , t * or one subexpression 
moving accross a power of two.) Note the value of b is 
effectly exact. For x outside <.9286194,1.857239)» a and b 
are non-zero and b dominates log x ? hence, the error bounds 
are 1 bit apart. For x In (.9286194,1,107239), log x 
collapses to 2t*p . But t=(y-c)/(y*c) where y-c is exact, y*c 
may lose half a bit, and the quotient involves further error. 
So those combine with the addition in 2t*p to make the total 
error. For x in (1.107239,1.857239), log x=((2t+p)+b with 
b=(1og c)/2 almost exactly, t and b may be of opposite sign. 


4, EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in the argument x, the error in the 
result is given approximately by e*/x. 


16 
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ROUTINE t AT&ti 


i. ROUTINE'S FUNCTION. 

1.1. Type. A FORTRAN external function. It accepts a floating¬ 
point arguBientf and returns a floating-point result. 

1.2. Purpose, To accept calls for ATAN by name at entry point ATAN 
and by value at entry point ATAN. • ATAN computes the inverse 
tangent function. 


2. METHOn. 


The input range to this routine is the collection of all definite 
in-range normalized floating-point quantities. The output range of 
this routine is included in the set of those floating-point 
quantities lying between -pi/2 and pi/? . 

The argument x is then transformed into an argument y in [0« 1/161 
by the range reduction formulae 

arctanCu) = -arctan(-u)» u negative* 

arctan(u) = pl/4 ♦ {pi/(i - arctan fl/u) I * u>l 

arctan(u) = arctanfk/16) ♦ arctan f(u - k/16)/(l ♦ u*k/16)l. 

where (l<u<l» and k is the greatest integer not 
exceeding 16*u. 

Finally arctanCyl (for y in 10. 1/161) is computed by the polynomial 
approximation s 

arctan(y) = y f a(l)*y3 a(2)*ys * a<3)*y'^ ♦ a(4)*y® 

where 

3(1) = -.33333333333312845, 
a(2) a .1999959958014464, 
a(3) a -.1428541305087450, 
a(4) a .11112281616126149. 

The coefficients of this polynomial are those of the minimax 
polynomial approximation of degree 3 to the function f over [0, 1/41 
where 

f(u2) = (arctan(u) - u)/u®. 

(The algorithm and constants are copyright 1970 by Krzysztof 
Frankowski, Computer Information and Control Science, University of 
Minnesota, 55455. Coding is by Larry Liddiard, University of 
Minnesota.) 


3. ERROR ANALYSIS. 

A graph of the relative error of approximation of the algorithm over 
to, 1/161 is shown in figure 7. The maximum absolute value of this 
relative error is 3.21)1 * 10**-16. An upper bound on the absolute 
value of relative error due to machine error has been established at 
4,761 * 10*»-13, Hence, an upper bound on the relative error in the 
routine is 4.764 ♦ 10»»-13 . 

For lOOO arguments chosen randomly from the following intervals, the 
following statistics on relative error were observed. 
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Interval 
from to 

Mean 

Standard 
Oeviation 

Mini mum 

Maxiffium 

-1. 

1. 

-1.58g£-17 

2.2i6E-15 

-6.823E-15 

5.539E-15 

-10. 

10. 

-2.348E-17 

1.940E-15 

-6.627F-15 

7.5fl5E-15 


4. EFPfCT OF ARGUMENT ERROR. 

If a small error e* occurs in the argument* the error in the result 
y is tiiven approximately by e*/fl ♦ y*l. 
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1. ROUTINE'S FUNCTION, 

1,1, Type, A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

!•’, Purpose, To accept calls by value for ATANH from FORTRAN 
programs, ATANH computes the inverse hyperbolic tangent, 

?, MFTHOO. 

The input range is the collection of all definite, in-range 
floating-point quantities in the interval (-1,0 , ♦•I • 01 • 

The range is reduced to til, II using the identity 

atanhf-x)=-at anh t xI, 

From the definition t anhf x ) = (eex-e^ f-x 1) / {eex ♦■e^ ^-x)) I one gets 
atanh(x) = 0,5»ln( (l+xl/fl-xit 

Usino the property I n (a ♦b 1= I n < a) t-1 n (h) , we can reduce the argument 
range of the above log to [,75,1,51 by extracting the appropriate 
mu! tip! e of In (?lJ 

at anh ( x) = 0. F'^n* I n ( ♦ 0 .5* I n ( ?♦ (-n I » (1 ♦ x) 7 (l-x> I 

Pewritinn the argument of log in the form (l+y)7(l-y), and 
substituting atanhfyi; 


?^(-n)»(l*x)-(l-x) 


atanhfx) = fl,5''‘n*-ln(?l*3tanh( 


-) 

2»(-n)»(l*x)*(l-x) 


This reduces th® range to I-0,2,♦n,?], 

The value of n such that f-n)*C1*xJ/(1-x) is in t,75,1,5) is the 
same as that such that ?♦ (-n) » f 1 *^xl / « 0.75* ( 1-x) I is in 11,2), If we 
write as a^2^i!i, a in [1,2), then 2»‘<-n-m) ♦ (1 + x)/a must be 
in [t,?>. If (i*x)>3 then -n'-m=0 and n=-m* If {l*x)<a then -n-m = l 
and n=l-fn. 

The function atanh (z> on t-0.2+0,2) is approximated by z+z^’p/q 
where p and a are 4th order even polynomials. The coefficients of p 
and g were derived from the (7th order odd)/(4th order even) minimax 
(relative error) ratioral form on t-0,2,+0,21 for atanh(z). 
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EPPOR ANALYSIS 


For abs(xl <0«2»n=0 and the torn is used and the error stays 

within the expected bound of 4.8E-1S. 

For abs<xJ>C,5, the term n*Cln<?)/2i dominates. This term is 
computed as n*C I n C2I/2-• 12*5)-n*.l 25-n*. 1?5 because the rounding 
error in representing ln{2)/2 is large? the above form makes the 
rounding error relatively small. Since n*.125 is exact and the 
dominating form, the two adds in (other) * 0 *.12'5«-n*. 125 dominate the 
error and the expected relative error of 8.3E-15 is the maximum 
observed error in this region. 

For 0.2<abs(x)<0.5,n=l and the term z={0.5*<1+x)-(1-x))/ 
( 0. 5** (1+x) (1-x) I may be relatively large. For abs(x)<0.25. the 

subtraction l-x=0.5-x+0.5 loses two bits of the original argument. 
In additiort z is negative in this range and some cancellation 
occurs in the final combination of terms. costing about one ulp. 
the actual upper bound in the region 0.2<abs(x)<0.25 is 19.4E-15, 
which is the overall upper bound. 

The errors aret 


rational form 
coefficient rounding 
round-of f 
upoer bound 
maximum observed 


fiECfill-iliLi® 

2.2 

< 0.1 

17.1 

19.4 

12.3 


4. EFFECT OF ARGUMENT ERROR. 


For small errors in the argument x* the amplification of absolute 
error is l/(l-x*l and that of relative error is x/C fl-x®) ♦atanh (x) ) . 
which increases from 1 at 0 and becomes arbitrarily targe near 1.0. 


e.g.. 18.8 

at 9.99 and 132 

at 

G.999, 

or 

approximate 1 

y -1/ 

( eps* 1 n (eps )) 

where x=l-eos. If 

X is 

known 

to 

more 

than 

single 

precision. th 

e following FORTRAN 

may 

be used 

to 

get a 

better 

result 


near l.Ol 
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OOUBLF X 


Ccompute XI 
SNGLX=X 

SHSNGLX=X-SNGLX 

Y=ATANH(SNGLXI ♦SHSNGLX/( (1+SNGLX) »SNGL (l-xin 


which is accurate 
accurate above 
ATANH(SNGLtXII, 


to sinqie precision 
this point, a 


for abs(x) and 

though stil! better 


1 ess 
than 


60498200 C 


23 





• 24 


60498200 C 


99990 POINTS. MEfiN R.E. 1.37E-IS RMS R.E. 2.S7E-15 










1 


»?OUTINE'S FUNCTION 


l.l. Type. A FORTRAN external function. It accepts an argument 
set comprising two floating-point arguments* and returns a 
floating-point result, 

1.?. Purpose, To accept calls for ftTAN2 by name at entry point 
ATANP and by value at entry point ATAN2. . ATAN2 computes the 
inverse tangent function of the ratio of two arguments. 


?. MFTHOn. 


The input range to this routine is the collection of all pairs (x*y) 
of definite in-range normalized floating-point quantities such that 
Cx.yl * (O.nt. 

The function ATAN?<xtV) is defined to be the angle (lying in (-pi* 
pin subtended at the origin by the point (yfX) and the first 
coordinate axis. 

The argument <x*y) is reduced to the first quadrant by the range 
reduct Ions 

ATAN2(x,y> = -AT AN2<-x ,v) , x<fi: 

ATAN2(x,y» = pi - ATaN2(x*-yl, x>0, y<0. 

The argument (x*y> is then reduced to the sector 
€(u,vS: u>0 & v<u & v>0> 
by the range reduction 

aTAN2(x,y) = pi/2 - ATAN2(y»xl, x>0 or y>0. 

Then ATANPCx.y) is evaluated as arctan(y/x)* using the algorithm 
described in +he method section of the routine ATAN "s description, 
(The algorithm and constants are copyright 1970 by Krzyztof 
Frankowski, Computer Information and Control Science* University of 
Minnesota* 55455. Codino is by Larry Liddiard* University of 
Minnesota.! 


T, ERROR ANALYSIS. 


See the error analysis of ATAN for properties of the algorithm used 
in computing arctan(y/x), 2000000 pairs of arguments (x,y) were 
randomly generated belonging to sets {:<u*v) s lul* Ivl < lO^'^kT, 
where k = -100, -99* ... * 100. The maximum absolute value of the 
relative error in the routine for these arguments was observed to be 
9 . 3^9 * 10**-15 for these random arguments. 

For 1000 arguments chosen randomly from the following intervals* the 
following statistics or relative error were observed. 
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In terval 
from 

O f X 

to 

Interval 

from 

of y 
to 

Mean 

Standard Minimum 

Deviation 

Maximum 

-1. 

1. 

-1. 

1. 

-3. 182E-rl6 

2.501E-15 

-l.GOlE-lt 

8. lElE-15 

-inn. 

lOO. 

-100, 

100. 

-2,4?qE-16 

2.E12E-15 

-1.012E-14 

8,374E-IF 


4, EFFECT OF flPGUHEWT ERROP. 


If small errors e(x> and e(yl occur in x and y respectively, the 
error in the result is given approximately by Cy*e(x) - x*e(yl)/(x2 
♦ y2) . 
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1. ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a complex 
argument and returns a fIoatlng-point result. 

1.?. Purpose. To accept calls by name and by value from FORTRAN 
programs for computation of the complex absolute value 
function. 


?. METHOD. 


The input range is the collection of all valid complex quantities 
whose absolute value does not exceed 1.?85*10®*® • 

Let X •• i»y be the argument. The algorithm used is* 

a. u •• maxdxltlyl}. 

V4>min(txl«lyll. 

b. If u or V fails a test for infinite or indefinltey go to 
step f. 

If u is zero* return zero to the calling program. 

c. r u/v 
w «- l*-r2 

t •- (T3/32 * 3/8Mw - 33/321 
= 3/8<ra * 87/32) 

<t is the initial linear approximation to <H-r*)**0.5) 

d. Heron's rule is applied in three stages, 
tfl) •• l/2(t ♦ H/t) 

t(2) l/2ltCl) ♦ w/t (D) 
t(3) «- l/2(t(2) ♦ w/t(2l) 

e. Return with u*t13) to the calling program if it Is not 
in finite. 

f. Call routine SYS=1ST to initiate error processing. 

g. Return to the calling program, unless a non-standard or 
fatal error recovery has been chosen for this routine. 

Note that a number of valid arguments are netted in step b. but 
these are returned to normal execution after further testing. 

Formulae used are 

|x»l*vl = SQRT(x+i*y» 

= maxClxl ,lyl)*Cl*r2)*».5 , 
where r = mint lx I.IyI)/max(Ixl.IyI). 

See the timing information in Appendix 0 for further details. 
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ERROR ANALYSIS 


The maxifflum absolute value of the error in approximating 
t<31 = SQRT(l+r2| by 

t = 33/32 ♦ 3/8ll+r2 - 33/32) 
tfl) = l/2(t ♦ (l*r2|/t) 
tC2l = l/2ftCl) ♦ (l+r*)/tCl)) 
t(3> = t/2<t(2) ♦ fl+r2)/tC2)) 

is 1.5306*10**-16» assumed when r=0« Hence an upper bound on the 
absolute value of error in the algorithm is 
1.53Qh*10**-16*max(lxl,|y1) « 

where x+iy is the argument. An upper bound on the absolute value of 
error in the routine due to machine round-off has been established 
at 8,512*10*^*-1*» * max<lxl«ly1) • Therefore, an upper bound on th 
absolute value of error in the routine is 8,527*10**14 
maxflxlflyl) , and an upper bound on the absolute value of relative 
error is 8.527 ♦ 10»*-14. 

For 10000 arguments chosen randomly from the interval 

the following statistics on relative error were observed, 

Nean Standard Hinlmum Maximum 

Oeviation 

-2.2qFE-lS 2,658E-15 -1,093E-1«» 5.967E-15 


4, EFFECT OF ARGUMENT ERROR. 


If a small error e(zl = e(x)*i*ety) occurs in the argument z = x*i*y 
• The error in the result u is given by 

e(u) = (xefx)♦yeCy))/u • 
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ROUTINE’S FUNCTION 


1.1. Type. A FORTRAN axTernaI function. It accepts 
arqument and returns a complex result. 

1.?. Purpose. To accept by name for CCOS from FORTRAN 
CCOS computes the complex cosine function. 


a complex 


programs. 


?. METHOD. 


If u and V are real numbers, then 

cos(u»-i.v) = cos (ul .coshivl-sin lu).sinhC v) .i . 

The argument is checked upon entry. The argument is invalid If 
either the real part or the imaginary part is infinite or 
indefinite, if the real part or the imaginary part is so large that 
precision will be lost during the computation, or if floating 
overflows occurs during the computation. If the argument is 
invalid, POS.INDEF. 4- i.POS.INCEF* is returned, and a diagnostic 
message is issued. If the argument is valid, C0S=5;iN is called at 
entry point COS,SIN for computation of the cosine and sine of the 
real part of the argument, and HYPER3. is called at entry point 
HYPERP. for computation of the hyperbolic cosine and sine of the 
imaginary part of the argument. The result is calculated according 
to the formula above and is returned to the calling program. 


T. ERROR ANALYSIS, 


The algorithm used in CCOS is the same as that used in CCOS. . See 
the description of CCOS. for the error analysis. 


U, EFFECT OF ARGUMENT ERROR, 


If a small argument error appears, then the error in the result is 
given approximately by multiplying the argument error by the 
negative of the complex sine of the argument. Hence, if a small 
error occurs in the complex argument and the error has absolute 
value e’, then the absolute value of the error in the result is 
given approximately by e' . (sinful® + sinhlv)®)♦*!/? , where u+i.v 
is the complex argument. 
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1. ROUTINE'S FUNCTION. 


1.1* Type. A FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by value for computation of the 
complex cosine. 


?. METHOD. 


The input range is the collection of all definite in-range complex 
quantities z = x ♦ l.y where lyl does rot exceed 741.67 and Ixl does 
not exceed • The formula used for computation is 

cosCz) = cosix + i.yl = costxl.cosh(y) - i.sin(x).sinh(y) . 
where x and y are floating-point Quantities. COS=SIN is called for 
computation of cosfxl and sin<x)« and HYPERR= is called for 
computation of coshfy) and sinhfy). The result is returned to the 
calling program - the real part in X6 and the imaginary part in X7 . 


7, ERROR ANALYSIS, 


(See the descriptions of COS=SIN anc HYPERB* for details.I If z « x 
* i.v is the argument* then the modulus of the error in the routine 
does not exceed 1.241 . 10»*l-13> ♦ 1.241 . 10**(-13l . expdyj). 

For 10900 arguments chosen randomly from the Interval t-l.*l.l*t- 
1..1.1. the following statistics on relative error were observed. 

Register Mean Standard Minimum Maximum 

Oeviat ion 

X6 -3.501E-15 3.827E-15 -1,413E-14 1.182E-14 

X7 -7.311E-15 g.884E-15 -5,059E-14 1.771E-14 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e(z) = e(xl ♦ i.efyl occurs in the argument z = x ♦ 
i.y* the error in the result is given approximately by 
-sin (z) .e (z} . 
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1. ROUTINE’S FUNCTION. 


1.1. Tyoe. A FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.2. Purpose. To accent calls by name for CEXP from FORTRAN 
programs. CEXP computes the complex exponential function. 


2. METHOD. 


If u and V are realt then 

expfu+i.w) = exp(u).cosfV) ♦ i.exofu).sin(v) . 

The argument is checked upon entry. It is invalid if the real part 
u or the imaginary part u is infinite or indefinite, if u is greater 
than 741.67 in absolute value* if v is so large as to lose accuracy 
during the calculation (i.e. v exceeds oi.2’>« in absolute value!* or 
if floating overflow occurs during the calculation. If the argument 
is invalid, POS.INOEF. ♦ i. POS.INOEF* is returned, and a 
diagnostic message is issued. If the argument is valid* the result 
is returned to the calling program. 


ERROR A4IALYSIS. 


The algorithm used in CEXP is the same as that used in CEXP. • See 
the description of CEXP, for the error analysis. 


4. EFFECT CF ARGUMENT ERROR. 


If a small error e* occurs in the argument u ♦ i.v* the error in the 
result is given approximately by e* • explu *■ i.v). Hence, the 
absolute value of the error in the result will be approximately 
le*|.exp(u). If the error in the argument is significant* the error 
in the result should be determined by substitution of possible 
argument values in the function. 
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ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by value for computation of tbe 
complex exponential function. 


?. HETHOO. 


The input range is the collection of all definite in—range complex 
quantities z = x ♦ i.y where |y| does not exceed pi.e'*® and |x| does 
not exceed 741.67. 

The formula used for computation is 

exp(z) = expfx + i.y) = expfx).cos(y) ♦ 1.exp(x).sin(y) 
where x and y are not floating-point quantities. 

COS=SIN is called for computation of cosCy) and sinty), and EXP. is 
called at entry point EXP. for computation of exp(x). The result 
is computed according to the formula and Is returned to the calling 
program. 


7. ER»OR ANALYSIS. 


(See the descriptions of COS=SIN and HYPEPB* for details.) If z = x 
♦i.y is the argument, then the modulus of the error in the routine 
does not exceed 1.378 . 10**(-13) » 1.378 . 10**(-13) . expdxl). 
If the real part of the argument is large, the error in the routine 
will be significant. 

For 10000 arguments chosen randomly from the following Interval, the 
following statistics on relative error of the components of the 
results were observed. 


Register Mean Standard Minimum Maximum 
to Deviation 

-1. 1. -1. 1. X6 

X7 

-670. 670. -2.210E14 2.2106E14 X6 

X7 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e(z) occurs in the argument z, the error in the 
result w is given approximately by w.e(z). 


-3.440E-15 3.784E-15 -1.428E-14 1.227E-14 
-5.831E-15 8.853E-15 -4.165E-14 1.242E-14 
-8.96ZE-15 4.669E-14 -3.176E-12 2.235E-14 
-1.071E-14 7.948E-14 -4.977E-12 3.723E-14 


Interval x Interval y 
from to from 
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ROUTINE t CLOG 


1. ROUTINE'S FUNCTION. 


1.1. Type, A FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by name for CLOG from FORTRAN 
programs. CLOG computes the complex logarithm function. 


p. METHOn. 


The argument is checked upon entry. The argument is invalid if the 
real or complex part is infinite or irdeflnitet or if both the real 
Dart and the complex part are zero. If the argument is invalid, a 
diagnostic message is written and POS.INOEF. ♦ i*POS.INOEF. is 
returned. Otherwise, CLOG, is called at entry point CLOG* for 
computation of the complex logarithm. The result is returned to the 
cal Iing program. 


T. ERROR ANALYSIS - see the description of CLOG* . 


4. Effect of argument error. 

If a small error e* occurs in the argument z, the error in the 
result is given approximately by e*/z. The modulus of this will 
give approximately the modulus of the error. 
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1. POUTINE*S FUNCTION, 

l.l* Type, a FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by value for computation of the 
complex logarithm, including converted calls from CLOG • 


?. METHOD. 


The input range to this routine is the collection of all definite 
in-range complex quantities which are non-zero, and whose absolute 
values do not exceed the largest floating-point number representable 
in the machine. 

The formula used to compute the complex logarithm is 
log z = logllzl) ♦ i.argizl, 

where |z| is the modulus of z. |z| is evaluated by routine CASS,, 
and the logarithm is evaluated by ALOG* • The function argCzl is 
evaluated by routine ATAN?, * argiz) always lies in the interval 
1-pi, pi) for z nonzero, definite and in-range. The result is 
returned to the calling program in XGeX? . 


3. ERROR ANALYSIS, 


Tests on a sample of lOflOHO random numbers distributed over the 
complex plane with distribution the product of two Cauchy 
distributions of zero mean returned a maximum absolute value for the 
relative error in the routine of 8.579 * 

For 10000 arguments chosen randomly from the interval 

the components of the results gave the following statistics 
on relative error. 


Register 

Mean 

Standard 
Devi at ion 

Minimum 

Maximum 

X6 

-7.1E0E-14 

4.603E-1E 

-4.435E-10 

4.213E-11 

X7 

-?.200E-16 

2.489E-15 

-1.114E-14 

8.085E-15 


4, EFFECT OF ARGUMENT ERROR. 


If a small error e(z) occurs in the argument z , the error in the 
result is given approximately by elzl/z. 
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, ROUTINE'S FUNCTION. 

1.1. Type. An auxilary routine of the math library. A floating¬ 
point argument is accepted and two floating-point results are 
returned. 

1.?. Purpose. To accept calls for COS.SIN from other routines* 
requiring the simultaneous computation of sine and cosine of 
the same argument. 


?. METHOD. 


Polynomials p(x) and a(x) of degrees 11 and 12 are used to compute 
sin(x) andcosfxl over the interval <-pi/4* pi/41* to which the 

argument is reduced. Upon entry, the argument x is multiplied by 2/ 
pi* and the nearest integer n to 2/pi . x is computed by double¬ 
precision addition of 2/pi . x to 200000000000000000003 * followed 

by rounded floating-point addition of the upper and lower halves of 
the result, n is normalized, and the argument x will cause a return 
of POS. INOEF. if the shift count in this normalization is zero' 
in other words, if x exceeds o i ,2'»« ( i . e., 221069929750 888.5...) in 
absolute value* Otherwise y = x n.pi/2 is computed in double- 
precision as the reduced argument for input to ply) and qly). 
sin(x) and cosix) are computed from these as indicated by the value 
mod(n*4). y lies in the interval <-pi/4*pi/4). The polynomials 
o(x) and q(x) are respectively 

s(n)x *■ slljxs ♦ s( 2 )xs ♦ sI3)xT ♦ sl4)x» ♦ s(5)x»» 

and 

c(0> ♦ C(l)x2 ♦ cC2)x‘» + cC3)x« ♦ cC4)x» ♦ c(5)xi« ♦ c(6)x** 
where the coefficients are given by 
SIQ) = .999999999999972 

s(i) = -.166666666685404 

s(2) = .833333331696029 . 10 •» -2 
s(7) = -.198412607353790 . 10 ** -3 
S(4I = .275548564509884 . 10 *♦ -5 
s(5) = -.247320720952463 . 10 *♦ -7 

cfOl = .999999999999996 

cll) = -.499999999999991 

cl 2) = .0416666666664705 

Cl3) = -.138888888698159 . 10 ♦♦ -2 

c(4) = .248015784673257 • 10 ♦♦ -4 
Ci5i = -.275552187277097 . 10 *» -6 
cl6) = .206291063476645 . 10 -8. 

These coefficients were obtained as follows. The polynomials of 
degrees 15 and 14 obtained by truncation of the MacLaurin series for 
sinix) and cos(x) were telescoped to form the polynomials pix) and 
qlx) of degrees 11 and 12. The method of telescoping polynomials 
fc.f., for example* C. Lanczos* Applied Analysis . 1956) consists of 
the (possibly repeated) removal of the leading term of polynomial by 
subtraction of an appropriate multiple of TIn)(a(X-x(0))) of the 
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same degree n, where 2/a is the length of the interval of 
approximation, and xCOJ is its centre. CTCnMxJ, the Chebyshew 
polynomial of degree n, is defined by TfnICx) = cos<n.arcos(x)I (1x1 
< 1) and satisfies the recurrence relation* 

T(OI(x) = 1 
T<1J fx» = X 

T(n+l){x) = 2xT(nMx) - T(n-ll(xl (n > II. 

T(n)(x> ffor n > 11 is the unique polynomial 2(n-l)*x**n ♦...of 

degree n whose maximum absolute valce over l-l, il is miniaial. This 
maximum absolute value is, of course, 1.1 

The formulae used for range reduction are* 
sinfxl = (-l)**n sin(y) 
cosixl s (-ll^^n coslyl 
if X = y ♦ n pi, n an integer* 
sin(xl“ cosfx - pi/2) 
cos(xl= -sinCx -pi/21 
if Pi/<»<x<oi/?. 

The input range is the collection of definite, in-range floating¬ 
point quantities whose absolute values do not exceed pi * E**® . 


ERROP ANALYSIS. 


The maximum absolute error in the approximation of sin(x} by p(x) 
over {-pi/*»,pi/4) is .1893 • 10 -14 and in the approximation of 
cos(x) by a(x) is .3687 • 10 ♦♦ -14 • Upper bounds on the machine 
round-off and truncation error over the input range (-pi/4,pi/4) 
have been established for p<x) at 7.523 • 10 -15 and for q(x) at 
1.401 . 10**-14 . Hence, the maximum absolute error and for qlx) a 
1.401 . 10 ♦♦ -14. Hence the maximum absolute error in this 
routine's computation of sine over (-pi/4,pi/4) is 9.416 • 10 * -15 
and of cosine is 1.770 ♦ 10 -14 • 


EFFECT OF ARGUMENT ERROR. 


Not applicable, since this routine is not directly called by the 
user's program. 
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1. ROUTINE-S FUNCTION. 


1.1. Type. fl FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by name for CSIN from FORTRAN 
programs. CSIN computes the sine function for complex 
arguments. 


2. METHOD. 


If X and y are real* than 

sirfx + i.yj = sin (x) .cosh(y) *■ i .cos (xl .sinhfy) . 

Upon entry* the argument is checked. It is invalid if the real part 
X or the imaginary part y is infinite or indefinite* if x or y is so 
large as to cause loss of precision in the calculation* or if 
floating overflow occurs during the calculation. If the argument is 
invalid, a diagnostic message is issued* and POS.INDEF. 
i.POS.INOEF. is returned. If the argument Is valid* the result of 
the computation is returned to the calling program. 


ERROR ANALYSIS. 


The algorithm used in CSIN is the same as that in CSIN, . See 
section T of the description of CSIN, for the error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


If a small argument error appears* then the error in the result is 
given aoproximately by multipling the argument error by the complex 
cosine of the argument. Hence* if a small error occurs in the 
complex argument and the error has absolute value e’* then the 
absolute value of the error in the result is given approximately by 
e* ♦ lcos(x)2 + sinh(y>?**l/2 * 

where x •• i*y is the complex argument. If the argument error is 
significant* the error in the result should be be found by 
substitution of the possible argument values in the function. 
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&Qyillg _l _CSIN. 


i. ROUTINE'S FUNCTION. 

1.1. Type. A FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by value for computation of the 
complex sine. 


METHOD. 


The Input range is the collection of all definite in-range complex 
quantities z = x + i.y where lyl does not exceed 741.67 and 1x1 does 
not exceed pi . . 

The formula used for computation is 

sinfx *■ i*y) = sinfxl * coshlyl ♦ i * cosfx) ♦ sinhfy) * 
where x and y are floating-point numbers. COS=SIN is called for 
computation of the cosine and sine of x, and HYPERB# is called for 
computation of the hyperbolic sine and cosine of y . The result is 
returned to the calling program - the real part in X6 and the 
imaginary part in X7 . 


7. ERROR ANALYSIS. 


(See the descriptions of HYPERB. and COSsSIN for details.) If z = x 
+ i.y is the argument, then the modulus of the error in the routine 
does not exceed 1.276 , 10»*(-13) ♦ 1.2R7 . 10»*(-13) . exp(|y|). 
The error in the routine is significant if the argument has a large 
(positive or negative) imaginary part. 

For 10(J00 arguments chosen randomly from the interval I-l.,l. ]*t- 
l.«l.], the following statistics on relative error were observed in 
the components of the results. 

Register Mean Standard Minimum Maximum 

Deviation 

X6 -5.592E-15 8.653E-15 -4.030E-14 1.228E-14 

X7 -4,97gE-15 5,877E-15 -3.165E-14 1,550E-14 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e(z) = e(x) * i.e(y} occurs in the argument z = x *■ 
i.y, the error in the result is given approximately by cos(z).e(z). 
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1. ROUTINf’S FUNCTION. 


1.1. Type, ^ FORTRAN external function. It accepts a complex 
argument and returns a complex result. 

1.?. Purpose. To accept calls by name for CSQRT from FORTRAN 
programs. CSQRT computes the complex square root function 
which maps to the right half of the complex plane. 


P, METHOD. 


For the algorithrat see the description of CSQRT* . Upon entry, the 
complex argument Is checked. The argument is invalid if its real 
part or its imaginary part is infinite or indefinite, or if floating 
overflow occurs during the calculation. If the argument is invalid, 
a diagnostic message is issued, and POS.INOEF. ♦ i^POS.INOEF. is 
returned. If the argument is valid, CSQRT, is called at entry point 
CSQRT. for the computation. The result is returned to the calling 
program. For the purposes of this computation, values returned by 
the routine will lie in the right half of the complex plane. 


T. ERROR ANALYSIS - see the description of CSQRT. . 


«♦. EFF«=^CT OF ARGUMENT ERROR. 


If a small error e* occurs in the argument z. The error in the 
result w is given approximately by e*/{2.w). The modulus of this 
will oive a approximately the modulus of the error. 
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1. ROUTINE-S FUNCTION. 

1.1. Type. A FORTRAN axternal function. It accepts a complex 
argument and returns a complex result. 

I.?. Purpose. To accept calls by value for CSQRT. from FORTRAN 
programs. CSQRT* computes the complex square root function. 


?. METHOD. 


The input range to this routine is the collection of all definite 
in-range non-zero complex quantities. If the argument is zero* zero 
is returned. 

II z X + i.y is the argument* then the result is given by w = u * 
i .V where u and v are determined as follows. 

Let a = (x2 4^ y2| ** 1/?, 

b = f(a *■ |x|l/2» »» 1/2* 

and 

c = y/<2.b). 

Then if x > 0* u = b and v = c* while if x <o* u = c . sigr(y) and v 
= b . s ign(y) . The result from this routine will always lie in 
the first or fourth quadrant of the complex plane* and complex 
quantities lying on the axis of negative reals will be taken by the 
routine to the axis of positive imaglnaries. 


ERROR ANALYSIS. 

The routine was tested against a sample of IQOOQQ random numbers 
distributed over the complex plane with distribution the product of 
two Cauchy distributions. The maximum observed modulus of relative 
error was 1,595 • IH^^f-lLI . 

For lOOPO arguments chosen randomly from the following interval, the 
following statistics on relative error of the components of the 
results were observed. 

Interval x Interval y Register Mean Standard Minimum Maximum 

from to from to Deviation 

-100. 100. -100. 100. X6 -4.79PE-16 2,652E-15 -9.774E-15 1.107E-14 

X7 -4.320E-16 2.655E-15 -9.726E-15 1.032E-14 

-lO.i®® 10.*®“ -10.*®® 10.*®® X6 -4.053E-19 2.632E-15 -1.012E-14 1.036E-14 

X7 -4,098E-16 2.637E-15 -9.520E-15 1.096E-14 
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EFPECT OF ARGUHENT ERROR, 


Tf a smal I error e(z) = e(xl i,e(y} occurs in the argument, 
error in the result w = u + i,v is given approximately by e(z)/ 
= Ce(x) ♦ i.ety)}/ 2(u ♦ i. v). 


the 
(2,w> 
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1. ROUTINE’S FUNCTION 


1.1. Type. FORTRAN exterral functions. It accepts a double¬ 
precision argument and returns a double-precision result. 

1.2. Purpose. To accept calls by value for computation of the 
inverse sine and inverse cosine functions. 


2. MFTHOn. 


The input range is the collection of all valid double precision 
guantities in the interval [-1.0.+1.01. Arguments outside this 
range will initiate error processing. 

The following identities are used to move the interval of 
approximation to t0,SQRT(.FI 1s 

arcsin(-x)=-arcsin(xl 

arccos(x)=pi/2-arcsin f x) 

arcsin(x)=arccos(sort(l-x*JI x>0 

arccos(x)=arcsin(sqrtCl-ya)) x >0 

Call the reduced value y. If y<.09375. no further reduction is 
performed. If not, the closest entry to y in a table of values 
^ arcs in < z) . s qr t f 1 -z 2 j , 2 = , 14 , , 39 , ,5 , 54 ) Is found, and the formula 

arcs in(x)=arcsin(z)+ arcsin(w) 

where w=x.sqrt(l-z^l-z.sqrt Cl-x^) 

is used. The value of w is in f-.07g2, + .08<i8) 

The arcsin of the reduced argument Is then found using a 15th order 
odd polynomial fwith quotient!* 

x+x3 Cc(3)+x2(c(5)+x2{c(7)♦x®Cc(11>+x®<c(13)+x®(c(151+a/(b-x®!)))!)! 

where all constants and arithmetic before cCll! are in double, and 
the rest is in single except the addition of c(ll), which has the 
form singIe+singIe=doub1e• The polynomial is derived from a minimax 
rational form (denominator is (b-x®)) for which the critical points 
have been perturbed slightly to make cdll fit in one word. 

To this value, arcsinfz! is added (from a table, and only if the 
last reduction above was done) and the sum is conditionally negated 
and 0,-pi/2, +pi/2, or pi is added to complete the unfolding. 
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ERROR ANALYSIS 


The maximum relative errros ares OASIN OACOS 

minlmax rational form error .f)A2E-29 .082E-29 

algorithm error 

(double precision coefficients) «76E-29 .48E-29 

maximum error observed 10.5E-29 12.5E-29 

The regions of worst error are (.09375».1446) for OASIN and 
f.989E,.9966) for OACOS. In these regions the final addition is of 
quantities of almost equal magnitude and opposite sign* and 
cancellation of about one bit occurs, the worst case being .1451- 
.9629. For OASIN. the polynomial range was extended to cover the 
region <•0821.•093751» where the worst error occurs. For OACOS. the 
extension is not used, so that the maximum relative error for either 
routine occurs in the region (.9956..9966) in OACOS. For 10000 
points randomly distributed in this region the maximum observed 
relative error in OACOS was 12.5E-29. 


4. EFFECT OF ARGUMENT ERROR. 


Tf a small error eps occurs in the argument x. the resulting errors 
in OASIN and OACOS are approximately eps7(1-x*)**.5 and -eps/(l- 
x2)**.5. The amplification of the relative error is approximately 
x7{T(x)»(l-xa)»*.5) where fix) is OASIN or OACOS. The error is 
attenuated for OASIN of abs(x)<0.75 and for OACOS of x>-.44» but may 
become serious for OASIN near -1 or +1 and OACOS near -1. If the 
argument is generated as l-y or y-1 then the identities 

asin(x)=acos(sqrt(l-x2)) 
acos(x)=asin(sqrt(l-x*)) 
asin(-x)=-asin(x) 
acos(-x)=pi*asin(x) 

may be used to get the full significance of y. When computing (1- 
x2) one should use a form such as (l-x*) = (H-x) ♦(l-x)=y*(2-y). 
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ROUTINE t^DAIAN. 

1, ROUTINE'S FUNCTION. 


1.1. Type, A FORTRAN external function. It accepts a double* 
precision argument and returns a double-precision result. 

1,?. Purpose. To accept calls from FTN compiled code for 
computation of the inverse tangent function. 


?. METHOn. 


The input range is the collection of all valid double-precision 
quantities. Other arguments will initiate error processing from 
OATAN= . Upon entry, the argument is loaded into registers XI and 
X? , and routine OATAN= is entered for ail remaining computations. 
See this routine's HETHOO description for further details. 


T, ERROR ANALYSIS. 


Sea section 6 of the description of routine OATAN* tor the error 
analysis. 


4. EFFECT OF ARGUMENT ERROR. 

See section 7 of the description of routine OATAN. 
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1. ROUTINE’ FUNCTION. 


l«l» Type. FORTRAN external function. The routine accepts a 
double-precision argument and returns a double-precision 
result. 

1.?. Purpose. To accept calls from FTN compiler code for 
computation of the inverse tangent function. 


?. method. 


The input range is the collection of all valid double-precision 
quantities. 

Computation is performed mainly in routine DATCOM. , and the 
constants used are listed there. 

a. Transfer return address from entry point word into B6 . 

b. Test first word of argument for infinite or indefinite. If 
either^ go to step i. 

c. 8t.-0. (B3 holds a mask HI.) 

87«-0« (B7 will hold closest multiple of pi/2 to absolute value 

of result.) 

d. R4 sign mask for argument. (84 holds MS , a mask for result’s 
sign.) 

e. XT x,t*- absolute value of argument. 

f. If absolute value of argument < 1. ♦ lump to routine DATCOM, at 
entry point DTN. to complete processing, 

g. X5 X3«- absolute value of argument, 

X4 Xl»l. 

BT .--0 

B7«-l 

h. Jump to routine DATCOM. at entry point DATCOM. to complete 
processing . 

i. Pick up parameter for error precessor. Call error processor, 
supplying given, argument and parameters. 

). If error processor returns control, return pi/2, with sign that 
stored in 64 . pi/2 is picked up by doubling an entry in a table 
starting at entry point ATN. in routine DATCOM. . 
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ERROR ANALYSIS 


10000 random arguments were generated in the interval 11/200.,?0Q.1, 
and the resulting graph of relative error versus argument is shown 
in the figure following this routine's description. In this sample, 
the- maximum absolute value of relative error is 7.183*10**C-29) . 
Groups of 40 double-precision arguments were chosen randomly in each 
of the following intervals, and the following statistics on relative 
error were observed. 


Interva I "s 
Lower 
Bound 

- 8 . 

.01 


Interval's 

Upper 

Bound 

8 . 

10 . 


Nean 

-1.995E-30 

-1.505E-30 


Standar d 
Deviation 

l.lOqE-29 
1.124E-29 


Minimum 

-2.063E-29 

-2.907E-29 


Maximum 


3.208E-2q 

2.745E-29 


The maximum absolute value of relative error in the algorithm is 
1.622E-29, and this occurs at 1.069781471095183. 


3.1. ALGORITHM ERROR. 


Up to 1/16, the plot shows the error in the economized 
polynomial* It is not centered because the first coefficient 
was forced to be 1. The Interval between (2n-l)/l6 and 
(2n+ll/16 is repeated twice (once reflectedi, but the waviness 
is damped because of adding atan(n/8) • The descrete }umps at 
(2n-ll/l6 are caused by the inaccuracies in atan(n/8) . Above 
1.0, the subranges are delimited by 16/(2n-l) • 

3.2. TOTAL ERROR 


Most of the errors can be traced back, with difficulty, to 
quirks in double precision addition. Note that the lower 
parts of the constants for pi and some of the atan(n/8l's are 
negative. While it allows the constant to be precise to an 
extra bit or two, the unpredictable sign wreaks havoc on the 
addition process. 


EFFECT OF ARGUMENT ERROR, 


If a small error e occurs in the argument x, the error in the result 
is given by e/IH-x^l. 
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1. ROUTINE'S FUNCTION, 


1.1. Type, ft FORTRAN external function. The routine accepts a 
doub!e-orecision argument and returns a double-precision 
result. 

1.?. Purpose, To accept calls from FTN compiled code for 
computation of the inverse tangent function with two 
arguments. 


?. METHOD. 


The input range is the collection of all pairs of valid double- 
precision auantities which are not both zero. Other arguments will 
initiate error processing from OflTAN?* , Upon entry, the arguments 
are loaded into registers XI, X?, X?, and X4 and routine OATAN?, is 
entered for all remaining computation. See this routine's METHOD 
description for further details. 


ERROR ANALYSIS. 


See section 6 of the description of routine 0ATAN2, for the error 
analysis* 


4. EFFECT OF ARGUMENT ERROR. 


See section 7 of the description of routine 0ATAN2* for the effect 
of argument error. 
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BQUTINE t OATAN?. 


1. ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function. The routine accepts an 
argument vector comprising two double-precision arguments, and 
returns a double-precision result. 

1.?. Purpose. To accept calls from FTN compiled code for 
computation of the inverse tangent function with two 
arguments. 


2. METHOD. 


The input domain is the collection of all pairs of valid double- 
precision quantities which are not both zero. 

Computation is performed mainly in routine OATCOM*. and the 
constants used are listed there. 

a. Test first words of both arguments to see if either is infinite 
or indefinite. If so. go to step }. 

b. Normalize first words of both arguments. 

c. Of first words of both arguments are zero, go to step i. 

d. 84 «■ sign mask of first word of first argument. 

R3 complement of sign mask of first word of second argument. 

88 •- return address in calling routine. 

87 - 1 . 

e. XF XT absolute value of first argument. 

X4 XI •• absolute value of second argument. 

f. X5 > X4 . Jump to routine OATCOM. at entry point OATCOM. to 
complete processing. 

g. XF <-> X4 
XT <-> XI 

Complement contents of 83. 

87 g, if first word of second argument is positive 
*- 2. if first word of second argument is negative. 

h. Jump to routine OATCOM. at entry point OATCOM. to complete 
processing. 

i. Supply message for "ARGUMENT VECTOR 0.0" . 

I . Pick up parameters for error processor. Cal I error processor, 
supplying given arguments and parameters. 
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k. If control returns from the error processor, return ♦INDEFINITE, 
to the calling orogram. 

3, ERROR ANALYSIS. 


A group of 40 random double-precision 
t.01,10.1 X t.01,10.]« and the following 
error were observed. 


arguments was chosen in 
statistics on relative 


Mean 

-?.64qE-30 


Standard 

Deviation 

2.161E-29 

* 


Min ifflum 
-6.108E-29 


Maximum 

3.115E-29 


The maximum absolute value of relative error in the algorithm is 
1.622F-29. 


4. EFFECT OF ARGUMENT ERROR. 

If small errors e* and e” occur in the arguments x and y 
respectively, the error in the result is given approximately by 

(X * e" - y ♦ e'J/fx® ♦ y2). 
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1. ROUTINE'S FUNCTION. 

1.1. Type. An auxifiary routine containing cofliiBon code for OATAN. 
and OATAN2* • 

1.?. Purpose. Coismon code for computation of OATAN and OATAN2 . 


P. METHOD. 


On 


entryt at both entry points OATCOM. and DTN. « 

3T = mask MI . 

Rk = mask MS = sign of final result. 

36 = return address after processing is complete. 

37 = closest multiple of pi/? to absolute value of result. 


In addition, at entry point OATCOM. , 

X4 XI = OU 

XE XT = OV , 

and at entry point OTN. , 

77 X3 = OU . 


Entry point ATN. is the start of an e ighteen-word table containing 
tan-i(n/8) (0Sn<8| in double-precision. Entry point OATCOM. 
corresponds to step a., and entry point DTN. corresponds to step 

b.. Constants used in the algorithm arei 

dT -.373 333 333 333 333 333 333 333 ?8E 915 
d5 .199 999 999 999 999 999 999 673 046 5?6 
d7 -.142 857 14? 857 14? 856 28r 180 055 289 
d9 .111 111 111 111 109 972 932 035 508 119 
cll = -.090 909 090 908 247 503 
Cl3 = .001 351 201 845 778 152 
a = -.085 666 743 757 593 039 
b = -1.133 579 709 202 919 6 

d3, d5, d7, d9 are double-precision constants, cll, cl3, a, b are 

single-precision constants. Arithmetic operations with d subscripts 
are done in double-precision, those with u subscripts are done in 
single-precision. Boolean operations have B subscripts. 

a. OQ *• OU/DV In double-precision. Carry OQ in X7 X3 . 

b. (OQ DA-OU at OTN. ) <Note that 0<no<l.) 

c. n *- nearest multiple of 1/8 to OQ . OL*-0 • 

d. If n=0 , go to step f. 

e. DA»-(OQ-N/8) / C1<-N/8»0A) , computed in double-precision. 
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f. If COA) (u)=0 t go to step h« 

XX «>(0AMu» *(u» COAKuJ 

X «-XX*»0.5 COAMul {(fnAMu) »(l) (DA) fuH(OA)} tu) 

♦{r} (OAXu))) 

O. DC XX *(01) «cl3 ♦Cdl XX »(dJ fd5 ♦((!» XX »(dl (d7 ♦Id) XX 

*<d) (d9 ♦Cd) XX ’(dl (dll ♦Id) XX »(uJ fcl3 +( 0 ) a/Cb -(u) 

XX)n ))H 

h, V - (OA)(u) ♦(dl DC *(d) ((DA)(u) -(d) (OA)(u) *(i) (OA)(u)/ 

((DA)(u) ♦(r) (OA)(u))) w v +(d) ((DA)(I) - X»((DA)(() ♦ 

{DA)(u) * (DA)(u)/ 

((DA) (u) Krl (OA) (u) ) J 

1* •- (87 * pi/2) -(8) B3 (upper and lower) 

1* £ *■ £ ♦(d) tan-*^(n/8)» tan-*(n/8) is obtained as a double- 

precision duantity from the look-up table* 

H, D «- (£ tld) w) -(B) (B3 -() BA) 

Xfi X7 «- p ♦ cleaned up. 

Return to address B6 by direct lump. 


3. ERROR ANALYSTS. 


Coefficients d3» d8, d7» d9. cll» cl3» at b were obtained by making 
the expression using these coefficients a minimax approximation to 
inverse tangent over 1-1/16,1/161, within the class of expressions 
obtained by varying these coefficients. (See descriptions of 
routines DATAN, and DATAN2, for error analyses.) 


A. EFFECT OF ARGUMENT ERROR. 


See descriptions of routines DATAN, and 0ATAN2, for effect of 
argument error. 
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POUTIHE'S FUNCTION 


1.1. Type. A FORTRAN external function. It accepts a double¬ 
precision argument and returns a double-precision result. 

1.2. Purpose. To accept calls by name for OCOS from FORTRAN 
programs. OCOS computes the cosine function. 


2. MFTHOh. 


See the description of OSNCOS* for the algorithm used in the 
computation. The argument is checked upon entry. It is invalid if 
infinite or indefinite or so large as to lose precision during the 
calculation. If the argument is invalid* POS. INOEF. is returned* 
and a diagnostic message is issued. If the argument is valid* 
OSNCOS* is called at entry point OCOS. for the computation. The 
result is returned to the calling program. 


EPPOR ANALYSIS - see the description of OSNCOS* • 


4. Effect of argument error - see the descriotion of OSNCOS* • 
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im t DCOSH 
ROUTINE’S FUNCTION 


1»1« Type. A FORTRAN external function. The routine accepts a 
double-precision argument* and returns a double-precision 
resuIt. 

1.?. Purpose. To accept calls from FTN compiled code for 
computation of the hyperbolic cosine function. 


?. METHOD. 


The input domain is the collection of all valid double-precision 
quantities whose absolute value is less than 1071*!og(2l. Arguments 
not in the domain will initiate error processing in routine DHYP» . 
Upon entry the argument is loaded into XI X2 before routine OHYP. is 
ailed. Csee the description of routine DHYF, for further details.) 


t. ERROR ANALYSTS. 


CSee the description of routine OHYP# for error analysis.) 
4. EFFECT OF ARGUMENT ERROR. 


fSee the description of routine OHYP# for effect of argument error.) 
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i, ROUTINE’S FUNCTION. 


i.l. 

Type. An 

auxiliary routine. 


1.2. 

Purpose. 

Common code for routines DEXP. * OHYP* « 

and 


DTANH . 




?. METHOD. 

Constants used in the routine ares 
1./loq(?} 

logf?) (in doubie-orecisionl 


d3 


.166 

666 

666 

666 

666 

666 

666 

666 

666 

709 

d5 

= 

. 833 

333 

333 

333 

33 3 

333 

333 

331 

234 

953E-2 

d7 


. 198 

412 

698 

412 

698 

412 

700 

466 

386 

658E-3 

dp 

= 

.275 

5 73 

192 

239 

858 

897 

408 

325 

90 8 

796E-5 

PC 

S 

-.47«< 

970 

880 

178 

988E-10 





pa 


.566 

228 

284 

957 

811E 

-7 






Db = 272.110 632 q03 710 
cll - .250 521 003 854 439E-7 


The algorithm is 

a. n •• nearest integer to x/fog 2 . 

y •- X - n ♦ log(2) . (Then y Is ir t-l/2*logC2), 172*1 og(2) 1 ) 

b. a ((y)(u) *(u) (y»(un**0.5 = (yMul (-(yj(u» *(1) 

q «■ (y)(u> ♦ (u) (y)(ul 

c. p q *(d) (d3 *(d) q *(dl (d5 ♦(d) q *(d) (d7 ♦(d) q *(d) 

(riq ♦(dJ q *(d) (cll ♦((!) q *(d) (pa / (pb-q» tpc) I) ) ) ) 


d. s - (y)(u) ♦(d) (y) (u) *(d) p 


e. (compute hm = SQRT{l^s**2) •) 
hi •- 3*0* ( (s) (ul ) 2 in single 
hi *- hi ♦ hi 
hk «• 2 * (l,+h!) 

hi «- ((y)(u) *(u) (y) (u) - h|)/hK-hi 

hm «■ b| ♦(u) (hk - (u) hi) *(u) (hl/hk) 

(hm now carries cosh-1, in sinale-precision) 


f. OS - s ♦(d) (((y)(l) ♦(r) (y)ll) *(u) hm) ttr) 

((s)(l) ♦(r) ((y)8u» *(l) (p)(u) ♦(r) (y) (u) *(r) (p)(l)))) 
(DS now contains sinh(y) in double-precision) 


g. DC hm ♦(d) (0S*ns-2*hm-hm*hm) / (2 (l.^hm) ) 

♦ ♦ 
evaluated in double 
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h. DX «- OS+OC 

i* Clean up DS» DC, g with 
X6 X7 •• DS 
Xf! XI *> DC 
X4 X5 •- OX 
n 

!• Direct Jump to 84 


T. EDPOR ANALYSIS. 


Not appticable. 

4. EFFECT CF ARGUMENT ERROR. 

Not applicable. 
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i. ROUTINE'S FUNCTION. 

l.i. Type. A FORTRAN external function. It accepts a double- 
precision argument and returns a double-precision result. 

1.?. Purpose. To accept calls from FTN compiled code for 
computation of the exponential function. 


2. METHOD. 


The input domain is the collection of all valid double-precision 
quantities lying in the interval 

t-975*log(2), 1070*log(2)l* i.e., t-675.84,741.671 . 

Arguments outside this range will initiate error processing from 
DEXP, . Upon entry, the argument is loaded into registers XI X2, 
and routine OEXP, is entered for the remaining computation. (See 
the description of routine OEXP# for further details.) 


3. ERROR ANALYSIS - see the description of OEXP, . 


4. EFFECT OF ARGUMENT ERROR - sea the description of OEXP, . 
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1. ROUTINE'S function, 

1.1. Type. A FORTRAN external function. It accepts a double- 
precision argument and returns a double-precision result. 

l.F. Purpose. To accept calls from FTN compiled code for 
computation of the exponential function. 


?. METHOn. 


The input domain is the collection of all valid double-precision 

ouantities lying in the interval f-975*log<2l.1070»logf2) • 

The argument reduction performed in routine OEULER, is 
X = <argument> 

V = X - n * log(2} 

where y = <reduced argument> is in t-1/2 log 2*1/2 log 21 
n is an integer. 

Most of the computation is performed in routine OEULER* « and the 

constants used are listed there. 

On input* the argument is in XI X2 * and on output, the result is in 

X6 X7 . 

a. Let X = <argument>. Save x in core. If 

l(xMuH > 177156400000000000008 * go to step g, 

b. Jump to routine OEULER* at entry point OEULER, with 34 = 
address for step c* X7 = upper part of x* X6 = lower part of x* 
XF = packed sign mask of x. 

On return from OEULER# , 83 = n * y4 = (OX)(u)* X5 = (0X)(l)* XO 
= OC)(ul, XI = 100(1), X6 = (00<u), X7 = (0S)(l), Here* n = 
nearest multiple of log 2 to x * y = x-n*log(2)* and OS 
Ssinhfyl* DC coshfy)-!, OX exp(y)-l* all in double-precision. 

c. w 1. +(0) (DC +(d) OS). Unpack w * increase exponents by n * 
and repack into X6 X7 . 

d. If upper word's exponent overflows, go to step g, 

e. If lower word's exponent underflows* go to step i. 

f. Return* with result in X6 X7 . 

g. Set parameters. Load up original argument. Call error 
processor. 

h. If error processor returns control* return. 
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i. Set parameters. Load up original argument. Call error 
processor• 

1. If error processor returns control return 0. 0. in X6 X7 . 

7. EPROR ANALYSIS. 


10000 random arguments were generated in the interval t-1/2 log S-* 
hf? log 23, and the resulting graph of relative error versus 
argument is shown in the figure following this routine's 
description. In this sample, the largest absolute value of relative 
error is 3.8S8E-29 . Groups of 100 double-precision arguments were 
chosen randomly in each of the following intervals, and the 
following statistics on relative error were observed. 


Interval•s 
Lower 
hound 
- 2 . 

-600. 


Interval"s 
Upper 
Bound 

?. 

700. 


Mean 

3.461E-31 

-8.631E-31 


Standard 
Deviation 

8.256E-30 

7.310E-30 


Minimum 

-2.632E-29 

-1.818E-29 


Maximum 

2.086E-29 

1.446E-29 


The approximation (described in the section on error analysis of 
routine OEULER, J is a mimimax approximation within the class 
obtained by varying the coefficients. 


3.1. ALGORITHM ERROR. 

The cunve for the algorithm error is barely distinguishable. 
It peaks at odd multiples of log 272 with a value of about 
.04E-29 . The algorithm error has essentially no effect on 
the total error. 

3.2. TOTAL ERROR. 

Except for fiddling with the exponent, OEXP ends with 1.0«-s 
where |s|<.3536 ? this addition is easy to do exactly when s 
is small and positive (see the plot lust above 0 and log 2). 
For s negative, the sum is less than 1, i.e., it crosses a 
band boundary, and it becomes difficult to produce an exact 
result (the plot shows exact or one bit low). When s< .25 
(e.g., .35<x<.45), it becomes even more difficult to prevent 
bits from dropping off in the low precision when lower sums 
overflow. 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e” occurs in the argument the error in the result y 
is given approximately by y*e' . 
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1. ROUTINE'S FUNCTION. 


1.1. Type. FORTRAN external functions. The routine accepts a 
double-precision argument at either entry point, and returns a 
double-precision result. 

1.?. Purpose. To accept calls from FTN compiled code for 
computation of the hyperbolic sine and cosine functions. 


?. HFTHOn. 


The input domain is the collection of all valid double-precision 
quantities lying in the interval (-1071*logC2)»1071*log(2)). 

Most of the computation is performed in routine DEULER. . and the 
constants used are listed there. The argument reduction performed 
in routine OEULER. is 
K = <argument> 
y = <reduced argument> 
y = x-r*log(?l 

where n is an integer, and y is in t-l/2*log(2),l/2*log(2)3 • The 
re-combination formula iss 
cosh(v+n*log 2) 

= (cosh(y)*sinh(yjI2**fn-1) ♦ (coshly)-sinhfy>)2**(-n-1) 
sinh(y+n*log 2) 

= (coshCy)♦sinh(y))2**(n-l) - (coshCy)-sinhly»I2»*(-n-1J 

At entry points OSINH. and OCOSH. , the argument is in XI X2 , and 
on exit, XF X7 holds the result. OSINH. corresponds to entry at 
step a., and OCOSH. corresponds to entry at step m. 

a. Let a = <argument> = XI X2 
b |a| Store b in X7 X6 . 

OF sign of a 

b. 8F packed zero 

B4 «- address of step g 


c. If (b>(u) < xmaxlu). Jump to routine OEULER. at entry point 

DEULER. • If (b)(ul > xmaxfu), go to step e. xmax is 
1071»log(2) • 

d. If (bl(l) < xmaxfl), lump to routine OEULER. at entry point 

OEULER. . 

e. XI X2 *■ a 

Set up parameters for error processor call with message 

"ARGUMENT TOO LARGE". If call was to entry point OCOSH. , 
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transfer contents of OCOSH* to DSTNH. • 

f. Call error processor. 

If (a)(ul is indefinite* return through entry point OSINH, with 
X7 = ♦INDEflMITE « Otherwise, return through OSINH, with S6 
X7 = ♦-INFINITE , the sign determined by 05 , 

g, C^?eturn from DEULER, with parameters 


03 = 

n 

X4 = 

(0X1(u) 

X5 = 

(Dxitn 

XO = 

(DC}(uJ 

XI = 

(DCKfl 

X5 = 

(DSI(uJ 

X7 = 

(OSKII 

where. 

if y = 1-n log(21, 

OX 

exp(y)-1 

DC 

cosh(v)-1 

OS = 

sinh(y) ,) 

If n=0 

, go to step 1 • 

If n>48 

, go to step k. 


u •- 2**(n-l) (OC+OSI in double 
V ♦ ?**(-n-ll tOC+OSI in double 
w 2**ln-l) ♦ u in double 
If n>24 , go to step h, 

w w $♦$ (2**<-n-llt^) (I > in double, sign $♦$ determined by 

85 , 

h« w *■ w <2**(-n-ll+vl lul in double, sign $♦$ determined by 

85 , 

i. X6 X7 w with sign same as sign of 05 • 

Return through entry point used to call routine, 

R. w ► tl,+OC*OS) ♦ 2*»fn-ll 
Go to step 1, 

l. If OSINH, entry, return through OSINH, • (Note that X6 X7 = OS 

,) 

X6 X7 «- 1, ♦ DC in double. 

Return through OCOSH, 

m. Let a *• XI S2 = <argument> 
b •- |a| Store b in X7 X6 , 

85 *- 1 

Go to step b. 


:7. ERROR ANALYSIS. 


10009 random arguments were generated in the Interval 
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t-1/2 fog 2,32 log 21 


for each of OSINH and OCOSH , and the resulting graphs of refatlwe 
error versus arguirent are shown in the figures following this 
routine’s description. In these samples, the maximum absolute 
values of relative error were 8.026E-29 for OSINH , and 4.4fl5E-29 
for ncoSH . The following statistics on relative error were 
observed in random samples of arguments in the intervals described. 


Entry 

Interval * s 

Interval 

•s Mean 

Standard 

Minimum 

Maximum 

Point 

Lower 

Upper 


Deviation 




Bound 

Bound 





OSINH. 

-2. 

2. 

8.516E-31 

1.086E-29 

-2.738E-29 

3.238E-29 


-600. 

700. 

-3.274E-31 

7.907E-3fl 

-2.645E-29 

1.651E-29 

OCOSH. 

-2. 

2. 

-2.055E-30 

1.217E-29 

-3.071E-29 

3.706E-29 


-600. 

700. 

-1.096E-3fl 

9.645E-30 

-2.733E-29 

1.904E-29 

3.1 . 

ALGORITHM 

ERROR. 






The curve 

for the 

algorithm error is barely distinguishable. 


Tt peaks 

at odd 

multiples of 

log 2/2 

with a value 

of about 


.04E-29 . The algorithm error has essentially no effect on 
the total error. 

OSINH 

The peaks are at odd multiples of log 2/2 below 33. • At 
A/.S'^log 2 , the algorithm error has a sudden peak* the reason 
is that it is at this point that the algorithm switches to 
OSINH(x)=expfx)/2 . This point was chosen because 2**{n-l) 
can be done correctly using an IX instruction to add n to the 
top of n.5. C48 would produce Indefinite.) Anyway, exp{x)/2 
is accurate enough. 

3.2. TOTAL ERROR. 

OCOSH 

The total error curves should be symmetric about the x=0 • 
The pattern shown should repeat until 47.S^log 2 (about 33.) 
at which point it will start looking like the OSINH and OCOSH 
curves. Between 8 and log 2/2 (=.3466), OCOSH is computed as 
1+c where 0<c<.75»SQRT(2)-1=.06066 t this is done fairly 
accurately, but the addition sometimes drops a bit in the low 
word. Above log 2/2 , the formula ends with a lot of addition 
and subtraction* for example OCOSSH (1.7443) = C 4 + 1/16)-4*.3*-sma I I 
stuff, where the ,3 is about what the sinh polynomial 
produced. Notice that the subtraction crosses a band and the 
exponent on 4».3 is only one less than the result; these facts 
make it difficult to keep from dropping bits. 

OSINH 

Do to log 2/2 , the error is predominated by the final add in 
the sinh polynomial. Just above log 2/2 the error is 


65 I 


60498200 C 




especially bad because of cancellation* Near log 2/2 , OSTNH 
is calculated via (1-1/41-s*-l/4*s where s is greater than 2'^*- 
2 and the result is less than 2**-l . The parts of the curve 
in the two ranges (.35,16.) and (16.,33.) have different 
shapes because of the shortcut taken in the latter range. 
(The split is at 23.5*tog 2 .1 Above 33. (47.5*log 2), the 
error curve is the same as for OEXP . 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in the argument x, the error in sinh(x) 
is approximately cosh{xl*e*, and the error in cosh(x) is 
approximately sinh(x)*e* . 
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1. ROUTINE’S FUNCTION, 


1.1. Type, A FORTRAN external function. It accepts a double¬ 
precision argument and returns a double-precision result, 

1.2. Purpose, To accept calls by name for DLOG from FORTRAN 
programs. DLOG computes the natural logarithm function. 


?, HETHOO, 


The algorithm used is given in the description of DLOG* , Upon 
entry, the argument is checked. The argument is Invalid if it is 
infinite or indefinitet or is not greater than zero. If the 
argument is infinite, indefinite or negative, POS, INOEF, is 
returned, Tf the argument is zero, NEG. INF, is returned. In any 
case, if the argument is invalid, a diagnostic message is issued. 
If the argument is valid, DLOG, is called at entry point DLOG. for 
the computation. The result Is returned to the calling program. 


E, ERROR ANALYSIS - see the description of DLOG, , 


•*, EFFECT OF ARGUMENT ERROR - see the description of DLOG. , 
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1. ROUTINE'S FUNCTION. 


i.l. Type. ft FORTRAN external function. It accepts a double- 
precision arg ument and returns a double-precision result. 

1.?. Purpose. To accept calls by value for OLNLOG. . calls 
generated by the use of OLOG or OLOGIO within FORTRAN 
programs. OLNLOG. computes the natural and common logarithm 
functions. 


?. HETHOn. 

The input range is the collection of all definite in-range double- 
precision ouantities which are greater than zero. 

Upon entry, the argument x is put into the form x = * w » where 

k is an integer, and 2**-l/2 S w < 2**l/2. Then log x is computed 
from 

log X = k . log 2 + log w • 

k • log2 is computed in double-precision, while log w is evaluated 
as follows. A polynomial approximation u is first evaluated in 
single-precision by 

u = cdJ.t + c(3).t3 + cfSI.ts ♦ c(7),t^ , 
t = (w - 11/(1 ♦ w) 

where the coefficients c(l), cl3), cl5) and cC7) are 
c(l) = 1.9R99gg993734000, 
c(3) = 0.6666694 SE.63 89k<», 
c(5) = 0.399657811051126, 
cC7) = 0.301005922238712. 

This approximates log with a relative error of absolute value at 
most 3,133 . 10**-8 over (2**-l/2 , 2*»l/2). Newton's rule for 
finding roots is then applied in two stages to the function exp(x) - 
w to yield the final approximation to log w. The two stages are 
aioebraically combined to yield the final approximation v* 

V = u - (l-x. exp (-u) I 

- Cl - X . exp <-u - (1 - X . exp(-ul))) . 

Writing z = 1 - x . exp(-u) , z is much less than 1, and v is 

computed by 

V = u - z(u) - z(l) - fz(u))2 . (.5 + zlu)/3) 

where z = z(u)*-z(() • This formula is obtained by neglecting terms 
which are not significant for double-precision. exp(-u) is 
evaluated in double-precision by the polynomial of degree 17 which 
is described in section 5 of the description of routine OEXP. . 

If entry was made at OLOGIO. , after k . log 2 *■ log w has been 

evaluated, the result is multiplied by log(lO) e in double¬ 
precision. The result is returned to the calling program. 


3, ERROR ANALYSIS. 

The maximum absolute value of the error of approximation of the 
algorithm to log x is 1.555 • 10**-29 over the interval [2**-l/2, 
2**1/21. A graph of the error in the algorithm versus argument is 
given in figure 16. An upper bound on the absolute value of the 
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machine round-off and truncation error (for arguments lying in 
1/?, has been established at 5.146 . 10*»-28. Hence the 

absolute value of the error in the routine over the interv al 
1/2. 2^»l/2] is not greater than 5.102 . 10»’>-28. The maximum 

absolute valut of the relative error of aoDroximation of the 
algorithm to log x over t2**-l/2. 2**1/21 is 2.266 . 10*»-28. An 
UDoer bound on the absolute value of the relative machine truncation 
and round-off error has been established at 1.486 . 10**-27. Hence 
an upoer bound on the absolute value of the relative error in the 
routine over the interval t2*»-l/2, 2*»1/21 is 1.713 . 10**-27. 


4. EFFFCT OF ARGUMENT ERROR. 

If a small error e* occurs in the argument x. the error in the 
result is given approximately by e*/x . 
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1. ROUTINE'S FUNCTION. 


1.1. Type. fl FORTRAN external function. It accepts a double- 
precision argument and returns a double-precision result. 

1.2. Purpose. To accept calls by name for OLOGIO from FORTRAN 
orograms. DLOGIO computes the common logarithm function. 


2. NETHOO, 


Upon entry, the argument is checked. It is invalid if it is 
infinite or indefinite, or if it is not greater than zero. If the 
argument is infinite or indefinite or negative, POS.INOEF, is 
returned. If the argument is zero, NEG. INF. is returned. In any 
case, if the argument is invalid, a diagnostic message is issued. 
If the argument is valid, OLNLOG= is called at entry point OLOGIO. 
for the computation. The result is returned to the calling program. 


T, ERROR ANALYSIS - see the description of DLNLOG . 


4, EFFECT OF ARGUNENT ERROR - see the description of DLNLOG . 
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ROUTINE X OMOO 
1. ROUTINE'S FUNCTION, 

1.1. Type. ft FORTRAN external function. It accepts two double¬ 
precision arguments and returns a double-precision result. 

1.?. Purpose. To accept calls by name from OMOD from FORTRAN 
programs. OMOD computes the modulus of an argument relative 
to a second argument. 


?. MFTHOO. 


The argument range is all valid double precision (x»y) such that [xT 
y]<?»96 and yjCQ. After argument checking, OMOO. is called to 
compute the result. The comparison tx/yl»2i'g6 is done by comparing 
exponents and, if necessary, coefficients. 


T. ERROR ANALYSIS - not applicable. 

4, EFFECT OF ARGUMENT ERROR - not applicable. 
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I SQUTINE_1 DMQD. 

1, ROUTINE'S FUNCTION 


Type. A FORTRAN external function. It accepts argument sets 
comprising two double-precision arguments* and returns double¬ 
precision results. 

1.?. Purpose. To accept calls by value for OHOO. * calls generated 
by the use of OHOD within FORTRAN programs. OHOO. computes 
the remainder of an argument relative to a second argument. 


2. METHOD. 


The argument range is all valid double-precision fx/yl such that tx/ 
yl<?^1070 and y70. The function computed by OHOO Cx*y) is 

x-tx/y]*y 

where ful denotes truncation. The value of x is repeatedly reduced 
by 4‘5-bit approximations to tx/yl until the reduced value lies in 
the range tE.signfy*x)). Since the result does not exceed 96 bits 
(see AMOD). the intermediate value of x does not exceed 98 bits and 
the reduction is done in triple precision* the result is always 
exact. 


7. ERROR ANALYSIS. 

Not applicable. Note that the only double-precision operations 
concerned in a determination of error are double-precision 
multiplication and double-precision subtraction. 4. EFFECT OF 
ARGUMENT ERROR - not applicable. 


74 


60498200 C 



,l-i2SI£L. 


i. ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a double- 
breclsion argument and returns a double-precision result. 

1.2. Purpose, To accept calls by name for OSIN from FORTRAN 
programs. OSIN computes the sine function. 


2. METHOO. 


The argument is checKed upon entry. It is invalid if it is infinite 
or indefinite or is so large as to lose accuracy during the 
computation. If the argument is invalid. POS. INOEF. is returned 
and a diagnostic message is issued. An argument Mill lose accuracy 
if it exceeds pi * 2^® in absolute value. If the argument is valid. 
OSNCOS. is called at entry point OSIN. for the computation. The 
result is returned to the calling program. 


I. ERROR ANALYSIS - see the description of OSNCOS. . 


V. EFFECT OF ARGUMENT ERROR - see the description of OSNCOS. • 
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1. !?OUTINE*S FUNCTION. 


1.1. Type. A FORTRAN external function. The routine accepts a 
double-precision arqutRent* and returns a double-precision 
result. 

1.?. Purpose. To accept calls from FTN compiled code for 
computation of the hyperbolic sine function. 


2, METHOD. 


The input range is the collection of all valid double-precision 
quantities whose absolute value is less than 107l*logf2). Arguments 
outside this range will initiate error processing in routine OHYP. . 
Upon entry, the argument is loaded into XI X2 , and routine OHYP, is 
called to complete the processing. (See the description of routine 
OHYP= for further details.) 


ERROR ANALYSIS, 


See the description of routine OHYP, for the error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


See the description of routine OHYP, for the effect of argument 
error. 
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1, ROUTINF’S FUNCTION. 

l.l. Type. FORTRAN external functions. The routine accepts a 
floubl e-precision argument and returns a doubi e-orecisicn 
resuf t. 

i.?. Purpose. To accept calls by value for DSNCOS. t calls 
generated by the use of OSIN or OCOS within FORTRAN programs. 
OSNCOS. computes the trigonometric sine and cosine functions. 


?. MFTHOO. 

The input range is the collection of all definite in-range double- 
precis i onquant it i es which are less than pi.?**® in absolute value. 
Upon entry* the argument x is made positive and is multiplied by 
?/oi in double-precision* and the nearest integer n to x • 2/pi is 
computed. At this stage* )x.2/pi| is checked to see that it does 
not exceed 2‘*‘^ . If it does* POS.INDEF. will be returned in X6 and 
a zero word in X7 • Otherwise* y = x - n . pi/2 is computed in 
double-precision as the reduced argumertt y lies in f-pi/4*pi/4). 

The value of niod(n*4)* the entry point called and the original sign 
of X determine whether a sine polynomial approximation pCx) or a 
cosine polynomial approximation qCx) is to be used* and also a flag 
to indicate the sign of the final result. 

The sine polynomial approximation is 

ofx) = a(l)x * a(3lxs * a(5)xs * a(7lx‘' + a<9)x» ♦ adDx^* t 
a(13)xi3 + a(15lx»^s ♦ adTIx^'? ♦ ad9)xi» + a(21)x2» 
and the cosine oolynomial approximation is 

a(x) = b(0) ♦ b(2)x2 ♦ b(4)x** ♦ b(6)x« ♦ b(8)x» ♦ b(10)x‘o f 
bd2)xiz ♦ bd4lx^‘» ♦ bCieixi® bd8)x»» + b{20)x2» 
for X in the interval (-pi/4* pi/4). 

The coefficients are 

ad) = .99999999999999999999999999999 

a(3) = -.16666666666666666666666666652 

a(5) = .83333333333333333333333270957 . 10 »» -2 

a(7) = -.19841269841269841269829134478 . 10 ♦» -3 

a(9) = .27557319223985890639440684401 . 10 »* -5 

adl) = -.25052108385441710113807647325 . 10 ** -7 

ad3) = .16059043836817941727119406461 • 10 ** -9 

ad5) = -.76471637307988608475534874891 . 10 -12 

3(17) = .281145706930018 . 10 »» -14 

a(19) = -.822042461317923 . 10 *» -17 

a(21l = .194362013130224 . 10 ♦♦ -19 

and 

b(0) = .99999999999999999999999999999 

b(2) = -.49999999999999999999999999919 

b(4) = .41666666666666666666666613902 

b(6) = -.13888888888888888888875543628 . 10*»(-2) 

b(0) = ,24801587301587301569992273730 . 10’^*(-4) 

b(10) = -.27557319223985877555866995711 . 10»»(-6) 

b(l?) = ,20876756987861921489874746135 . 10**C-8) 

bC14) = -.11470745595858431549595076575 . 10**(-10) 
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bU6) = .4779476968223*3311593310626721 . 10»*(-13) 
bUS) = -.156187668345316 . 10**(-15> 

b(20) = .408023947777860 . 10**C-18) . 

These polynomials are evaluated from right to left in double¬ 
precision using an in-stack loop. The sign flag is used to give the 
result the correct sign, before return to the calling program. 


3. FRPOR ANALYSIS. 

Graphs of the errors in approximating sinfxl and cosfxl by p(x) and 
qCx) over the interval C-Di/4,pi/4) are given in figures 13 and 14. 
The maximum absolute value of the error of approximation of o<x) to 
sinixl over <-pi/4,pi/4) is ,2570 . 10»*<-28), and of q(x» to cosCxI 
is ,^786 . 10**(-28), Upper bounds on the machine round-off and 
truncation error over (-pi 74,pi/4) have been established for p(x) at 
l.'^kT . 10*^*(-27) and for q(x) at 1.364 . 10»*(-27) . Hence an 
upper bound for the absolute value of error on this routine’s 
computation of sine over (-pi/4 , pi/4) is 1.769 . 10*»(-27) and of 
cosine is 1.402 . 10**(-27). 


4. EFFECT OF APGUMENT ERROR. 

If a small error e* occurs in the argument x, the resulting error in 
sin is given approximately by e" , cos(x) , The resulting error in 
cos is given approximately by -e- . sin(x). If the error e* becomes 
significant, the addition formulae for sin and cos should be used to 
compute the error in the result. 
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1. ROUTINE’S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a double- 
precision argument and returns a double-precision result. 

1.?. Purpose. To accept calls by name for OSQRT from FORTRAN 
programs. DSQRT computes the square-root function. 


?. MFTHOn. 


The argument is checked upon entry. It is invalid if it is 
infinite, indefinite or negative. If the argument is invalid, POS. 
INOFF, is returned, and a diagnostic message is issued. Otherwise, 
nS0PT= is called at entry point OSQRT. for the computation. The 
result is returned to the calling program. 


error ANALYSIS - see the description of OSQRT. , 


4. effect of argument error - see the description of DSQRT. . 
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1. ROUTIME'S FUNCTION. 

l.l. Type. A FORTRAN external function. It accepts a double- 
precision argument* and returns a double-precision result. 

i.F. Purpose. To accept calls by value for OSQRT. « calls 
generated by the use of OFQRT within FORTRAN programs. OSQRT# 
computes the square root function. 


?. MFTHOO. 

The input range is the collection of all double-precision quantities 
which are ^ero or positive* and are in-range and definite. Upon 
entry, the argument x is checked for a zero value. If it is zero, 
zero is returned. Otherwise the argument is put into the form 
X = 2 (2.el . y, 

where e is an integer, and .2F<y<l,Q0. The result returned is 2**^6 
* and y**(l/2) is calculated as follows. A fourth-order 

Chabyshev approximation pCy) is evaluated to obtain a single- 
precision initial approximation to , where y is the upper 

half of the double-precision argument. Heron’s rule ( z (n*!) = (z (n) *■ 
V/7(nl)/2) is applied in two stages in single-precision to give a 
single-precision approximation, and this is followed by an 
application of Heron’s rule in double-precision to give the final 
double-precision approximation. 

The polynomial pCyl is 

o(y) = . 182481f»34fl43495 

1.5462934655996 . y 

- 1.4758658070997 . y* 

T 1.06285652589999 . y^ 

- .323987345020001 . v^. 

This is evaluated as 

p(y) = cl((y T s* ^ + ly «■ S> iillly ♦ S)® ♦ £l 

where 

£ = -1.070137377460206 
51 = -1.391599471253464 
= 2. 286166868052419 
£ = -.00601995587532198 
c = -.323987345020001 

The final result is packed with the correct exponent e, and returned 
to the calling program. 


3. ERROR ANALYSIS. 


The algorithm error is at most 2.05E-31 falways positive). The 
round-off error in computing the single-precision approximation x is 
exactly 1/2 ulp. 

Including algorithm error, x may have Just over 1/2 ulp error, so x^ 
may have Just over 1 ulp error? since x is an approximation to the 
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single part onty« the total error in x® may exceed 2 ulp when 
y>x^lmax» 7.55E-15I* Then y-x* may contain 50 significant bits* and 
the error range for y-x* is (-1.78E-15,+3.55E-15)* and that for (y- 
is (-8•88E”15*+3»55E-15I. Relative to x* this error is C- 
&•71E-29*♦?.68E-291• In order to get this error* the error in x 
must be at least 7,11E-15* so the resulting error after the last 
Heron step is in <2.52E-2g,2.85E-29). The total error is in C- 
4»18E-29* 4-5.55E-29I . The maximum observed error in 100000 points 
randomly chosen in [1*4) was 3.19E-29'? the maximum in 200000 points 
randomly chosen in Cl.0,1.5) was 3.89E-29. 


4. EFFECT OF URGUMENT ERROR. 

If a small error e* occurs in the argument x* the error in the 
result y is given approximately by e*/(2.y). 
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PQUTINE t DTAN 
1. ROUTINE'S FUNCTION 


1»1. Type* A FORTRAN external function* It accepts a double- 
precision argument and returns a doubie-orecision result. 

1.2. Purpose. To accept calls by name for OTAN from FORTRAN 
programs. OTAN computes the trignonmetric tangent function. 


?. NETHOO. 


The argument is checked upon entry. It is invalid if it is 
infinite, indefinite or negative. If the argument Is invalid* 
ROS. TNOEF. is returned* and a diagnostic message is issued. 
Otherwise* OTAN= is catted at entry point OTAN. for the computation. 
The result is returned to the calling program* 


ERROR ANALYSIS - see the description of DTAN* 


k. effect of ARGUMENT ERROR - see the description of OTAN. • 
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SIMLl£j£_iJnAiU 

1. ROUTINE’S FUNCTION 


i.l. 


Type. A FORTRAN 
precision argument. 


external functiop. It accepts a double- 
and re'turns a doubi e-orec ision result. 


1.?. Purpose. To accept calls by value for OTAN* , calls generated 
by the use of OTAN within FORTRAN programs. OTAN. computes 
the trigonometric tangent function. 


?. METHOD. 


The input range is the collection of all double-precision quantities 
which are zero or positive, and are in-range and definite. Upon 
entry, the argument x is checked for a zero value. If it is zero, 
zero is returned. Otherwise the argument is put into the form 
X = 2 ** C2.eJ . y. 


where e is an integer, and .25?y?1.0G. The result returned is 
* Cy)**Cl/2l, and y**(l/2) is calculated as follows. A fourth-order 
Chebyshev approximation o(y) is evaluated to obtain a single- 
precision initial approximation to y**l/2, where y is the upper half 
of the double-precision argument. Heron’s rule (zCn+1)=(z(n) + y/ 
z(n>l/2» is applied in two stages in single-precision to give a 
single-precision approximation, and this is followed by an 
application of Heron’s rule in double-precision to give the final 
double-precision approximation. 

The polynomial pfyl is 

oCyJ * .182481834043‘*gF 

* 1.5462934655996 . y 

- 1,4758658070997 . yZ 
’ 1.06285652589999 . y3 

- .323987345020001 . y**. 

This is evaluated as 

p(yl = c(((y ♦ e * (y * el* b ><fy *• el®*- a)* c) 

where 


e = -1.070137377460206 
a = -1.391599471753464 
b = 2.286166868052419 
c = -.00601995587532198 
c = -.323987345020001 

The final result is packed with the correct exponent e, and returned 
to the calling program. 

A graph of the relative error in the algorithm of approximation to 
square root in doubie-orecisior over 1.25, 1.001 is given in figure 
12. The maximum absolute value of the relative error of 
approximation of the algorithm is 1.230, 10**(-31)» An upper bound 
for the relative error due to machire round-off and truncation 
errors has been established at 2.524 . 10*»(-28) . Hence the 
absolute value of the relative error ir the routine is less than or 
equal to 2.524 . 10*»{-28» . 
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3, ERROR ANALYSIS. 
h, EFFECT OF ARGUMENT ERROR. 

If a small error e* occurs in the argument x. the error in the 
result y is given approximately by e^^'f^.y). 
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1. WOUTINP-S FUNCTION. 


1.1» Type* A FORTRAN external function. The routine accepts a 
double-precision argument^ and returns a double-precision 
resut t, 

!•?. Purpose. To accept calls from FTN compiled code for 
computation of the hyperbolic tangent function. 


2. METHOD. 


The input 
quant ities* 
process ing 
into XI X? 
computation 
detaiIs.l 


domain is the collection of all valid double-precision 
Arguments outside the domain will initiate error 
in routine OTANH. . Upon entryt the argument is loaded 
f and routine OTANH. is entered to complete the 
(See the description of routine OTANH* for further 


T. ERPOR ANALYSIS. 


See the description of routine OTANH. for the error analysis. 

4. EFFECT OF ARGUMENT ERROR. 

See the description of routine OTANH. for the effect of argument 
error. 
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1. ROUTINE'S FUNCTION, 


1.1. Type. A FORTRAN external function. The routine accepts a 
double-precision argument, and returns a double-precision 
result, 

1.2. Purpose. To accept calls from FTN compiled code for 
computation of the hyperbolic tangent function. 


?. NFTHOO. 


The input domain is the collection of all valid double-precision 
quantities. Arguments outside the domain which are indefinite will 
initiate error processing. Most of the computation is performed in 
routine OEULFR, , and the constants used are listed there. The 
argument reduction performed is* 

<il argument in (-47*log 2,47*log 2) but not in 
(-17*log 2,l/2*!og 2) 

X = <argument> 

V = <reduced argufflent> 
y = 2x - n ♦ log 2 

where n is an integer, and y is in f-l/*loq 2,l/2*Iog 2J 
tanh(x) = u/v where 

u s 1 - 2»*-n - 2**-n * (OC-OS) 

V = 1 - 2**-n ♦ 2**-n * <OC-OS) 

(iiJ argument in (-l/2»log 2,lz»log 2) 

X = <argument> 

V = <reduced argument> 
y = X 

tanh(x> = 0S7I2.*0C) 

(ili) argument outside <-47»log 2,47*log 2) 

X = <argument> 
y = <reduced argument> 

tanhCxI = 1 - 2((1+OC-OSI * 2**-n - (d + DC-DS) * 2**-n)2) 

In Ci), (ii), and Ciii), OC cosh<y)-l and OS slnh(y). 

On entry to OTANH. , 1 X2 holds the argument, and on exit, X6 X7 
holds the result. 

a. Let a = XI X2 = <argument> 

X7 X6 •- b «■ |a| 

RE •- sign mash of a 
XR •- packed zero 

- 1 

R4 •- address of step e 

If exponent of first word of a is <-4g, direct Jump to routine 
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OEULER. at entry point OEULER. . 

X7 X7 ♦ 2 
X6 XE * 2 

84 •- address of step c 

If exponent of first word of a is <-42, direct Jump to routine 
DEULER* at entry point OEULER* • 

b. X6 X7 *■ $<-$1* with sign obtained from 85 
If a is definite* return* 

Set parameters for a call to error processor* 

Call error processor* 

If control returns from error processor, return* 

c. fOn return from OEULER* , 

= n = integer offset in argument reduction, 

X7 X6 = n*log2 ♦ y 
X4 = COXI<u) 

X5 = fOX)(1) 

XO a fOCJ fu) 

XI a fOCMI) 

X5 a fOSJ <ul 
X7 a <nSIIII 
where 

OX exofyl-l 
DC coshfyl-l 
ns sinh(y) J 
Tf n > 47 , go to step f* 
u l*-2**-n - 2**-n (OC-OS) 

V •- l.*2**-n ♦ 2**-n <DC-nS) 

d* we u/v , in double 
Go to step g* 

e* u e OS 

V e i.eoC , in double 
Go to step d* 

f* w e 1, - 2 * <<1*+0C-DSI * 2**-n - (C1.*0C-DS> ♦ 2*»-n>2» 

fevaluated in double, although only second word of 1* is 
affected! 

g* Clean up w , affix sign in 85 , and leave in X6 X7 * 

Return* 


ERROR ANALYSIS. 


lOOOO random arguments were generated in the interval 
t-l/2»log 2,3/2*log 21, 

and the resulting graph of relative error versus argument is shown 
in the figure following this routine*s description* In this sample, 
the maximum absolute value of the relative error is 8*58iE-29 * 
Random samples of 100 arguments were generated in the intervals 
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1isted, 
observed 


relative error 


were 


and the following statistics on 


Interval's 

Interval•s 

Mean 

Standard 

Minimum 

Maximum 

Lower 

Upper 


Deviation 



Bound 

Bound 





-2. 

2. 

3.O11E-30 

1.7353-29 

-6,675E-29 

7.436E-29 

-30. 

30. 

1.640E-30 

9.7603-30 

-3.692E-2g 

2.544E-29 


3.1. ALGOPITHM ERROR 


The algorithra error is insignificant. It is predominated by 
the error in the sinh expression in OEULER, . but by various 
folding actions, the error is damped even further. 

3.2 TOTAL ERROR 


The error plot should be symmetric about the origin. In the 
range (n,.5l the error is dominated by the code to divide s/ 
secondarily are the errors in s and in adding l*-c . 
Just above .5 . several factors conspire to create errors: an 
addition of numbers of opposite sign in the numerator, an 
addition in the denominator, and a division. (The errors in 
evaluating sinh are partially damped and fairly insignificant 
in comparison.) Up to 16.5 to 16.5 C23.75*log 21. the result 
is slightly less than l.O and the error is almost totally due 
to imprecise division of slightly imprecise arguments. From 
16e5 to 64.0 (2®). the result is perfect because it is 1- 
(iower precision stuffi where the was computed in double¬ 
precision. Above 64.0 (not shown), the error will taper off 
to zero because the arswer will be 1.0 while the true value is 
closer to 1.0 than 2**-96 . 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in 
result is oiven aoproxlmately by e* 


the argument x. 
* sech2(x) . 


the error in the 
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ROUTINE-S FUNCTION. 


1.1. Type. ft FORTRAN exponentiating routine. It accepts an 
argument set comprising two double-precision arguments, and 
returns a double-precision result. 

1.?. Purpose, To accept calls by name for OTOO* , generated by 
FORTRAN programs which raise double-precision ouantities to 
double-precision exponents. 


?. HFTHOO. 


The result is calculated byt 

result = expl exponent . iog(base)l. 

Upon entry, the argument set is checked. It is invalid if either 
argument is infinite or indefinite, if the base is negative, or if 
the base is zero and the exponent is not greater than zero, or if 
floating overflow occurs during the computation. If the argument 
set is invalid, POS.INOEF. is returned, and a diagnostic message is 
issued. Otherwise, DTOO* computes the result according to the 
equation above. The result is returned to the calling program. 


3. ERROR ANALYSIS. 


The algorithm used in routine DTOD* is the same as that used in 
routine OTOO. , the cal I-by-value counterpart. See section 3 of the 
description of OTOO. for the error analysis. 


4. effect of argument ERROR. 


If a small error e" occurs In the base b and a small error e*" 
occurs in the exponent p, the error in the result is given 
approximately by 

b D . (p/b.e* ♦ log(bl.e*"». 

The absolute error is approximately the absolute value of this 
expression. If the errors in the arguments are significant, the 
error in the result should be found by substitution of the possible 
argument values in the expression b ** p. 
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1. ROUTINF'S FUNCTION. 

1.1. Type. A FORTRAN exponentiating routine. It accepts an 
argument set comprising two double-precision arguments^ and 
returns a double-precision result. 

1.2. Purpose. To accept calls by value for DTOD, « calls generated 
by FORTRAN programs wMcb raise double-precision bases to 
double-precision exponents. 


2. MFTHOO. 

The input range is the collection of all argument sets (b.p) for 
which b and p are definite in-range double-precision quantities, b 
is positive, if b is zero then p is greater than zero, and b**p is 
in-range. 

The formula used isi 

b*»p = exo(p * log bl 

where b > n. Upon entry, DLNtOG. computes Jog b, and DEXP, computes 
expfp.tog bl . The result is returned to the calling program. 


t. ERROR ANALYSIS. 

10,0fltl pairs of double-precision random numbers were generated, with 
distribution the product of uniform distributions over (.5, 1.5) and 
(-in, 10). The error in the routine's computation of b**p was 
determined for each of these pairs. The maximum absolute value of 
the relative error in this routine for these 10,000 pairs was found 
to be 2.977 » 10*»{-25). 


4. EFFFCT OF ARGUMENT ERROR. 

If a small error eCbl occurs in the base b and a small error e(p) 
occurs in the exponent p, the error in the result r is given 
approximately by* 

r ♦ (log b » elp) + p * e(b)/b) . 
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1. ROUTINE’S FUNCTION. 


1.1. Type. A FORTRAN exponentlat in<| routine. It accepts an 
argument set consisting of a double-precision argument and a 
fixed-point arqument. and returns a double-precision result. 

1.?. Purpose. To accept calls by name for OTOI* from FORTRAN 
programs, generated when the programs raise double-precision 
quantities to fixed-point exponents. 


?. METHOD. 


The argument set is checked upon entry. It is invalid If either 
argument is infinite or indefinite* or If the base is zero and the 
exoonert is not greater than zero. If the argument set is invalid* 
a diagnostic message is issued* and POS. INDEF. is returned. 
Otherwise, DTOI* is called at entry point DTOI. for the 
computation. The result is returned to the calling program. 


ERROR analysis - not applicable. 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in the base b* the error in the result 
will be given approximately by r * b**(n-l) • e’, where n is the 
exponent given to the routine. 
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1. ROUTINE'S FUNCTION. 


1.1. Type, A FORTRAN exponentiating routine. It accepts an 
argument set comprising a double-precision base and a fixed- 
point exponent, and returns a double-precision result, 

1,?. Purpose, To accept calls by value for OTOI* , generated by 
FORTRAN programs which raise double-precision guantities to 
fixed-point exponents. 


?, METHOD, 


Let b be the base and p f?0J the exponent. If p has binary 
representation 0 00 . . . 01 Cn ) i <n-l) . . . i 11 U (0) where each il))C0<|<n) 
is 0 or 1, then 

D = i(0».2o ♦ i(l).2‘ *...♦ ifnl,2 ** n 
and n = tlog(2)p] = greatest Integer not exceeding log(2)p. Then 
b »* p = Prod €b 2 *♦ ) s 0 < 1 < r 1 iIJ) = 1> , 

The numbers b = »» 2 »» o , b^, b^,..,, b ** 2 n are generated by 
successive sauarings, and the coefficients i10), ,..,i(n) are 

obtained as the sign bits of successive circular right shifts of o 
within the computer. A running product is formed during the 

computation, so that smaller powers of b and earlier coefficients 
i(1 I may be discarded. Thus, the computation becomes an iteration 
of the algorithm 

b**D=lifp=fl 

b o = lb®) p/2 if p > Q and p is even 

b p = b . Ib2) ** Cp-lJ/2 if p > 0 and p is odd. 

Upon entry, if the exponent p is negative, p is replaced by -p and b 

is replaced by 1/b. b is double-precision, say b = x<u)*x<l). 1/b 

= <1/bl(u)»f1/bl11) is given in terms of x(u) and x(l) by the 

formulae below, where n is the normalization operation and the 
subscript I on one of the operations +, -, and . indicates that the 
coefficient of the result is taken from the lower 48 bits of the 96 
bit result register, and the exponent is 48 less than the single¬ 
precision coeffient's exponent, 

(l/b)(u) = n(l/x(u) + <(fnll-fl/xfu») .xCu) 

II -III Cl/x<u).x(u)) - Cl/x(u) .III x(ul) 

- Il^xlu) . x(l}/xlu))l 

♦ (l/x(ul ♦III lllnf l-(l/xlul) . x(un 

+ Cl -III ll/xlul) , xluin - Cl/x(u)) .III x(u» 

- (l/x(ul) . xlin/xlul) 

and 

(1/blIII = nl...) +111 <...! . 

In the routine, doubte-precis ion' quant ities x = xlu)*xll) and y = 
ylul^yll) are multiplied according to 
x.y = lx.yl (u)♦CX.yl(I ) 

where 

Ix.yKul = ( I fxlu) .y( f 1} + fxlll.yluM) 

♦ (xCu) .III ylulIJ ♦ IxIul.yCuM . 
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and 

(x.ytd) = f(Cx(ul . yCit) * fx(tl . yCuM) 

» (x(u) .fll yfu)l) (x(u}.y(un . 

The Input range is the col lection of pairs of arguments b» p for 
which p>0 if b is zero, all quantities are definite and in-range, 
and the result is ir-range« 


T. ERROR ANALYSIS - not applicable. 


%. EFFECT OF ARGUMENT ERROR. 


Tf a small error e* occurs in the base b, then the error In the 
result is given approximately by p • b**lp-l) * e*, where p is the 
exponent. Tf the error e* is significant, the absolute error in the 
result is bounded above by 

Ipl * maxdbl , lb ♦ e*ll**Cp-ll * le*l. 
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l, ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN exponentiating routine. It accepts an 

argument set consisting of a double-precision argument and a 
floating point argument, and returns a double-precision 
r esu It. 

1.2. Purpose. To accept calls by name for OTOX* generated by 

FORTRAN programs raising double-precision quantities to 

floating- point exponents. 


?. METHOD. 


The argument set is checked upon entry. It is invalid if either 
argument is Infinite or indefinite, if the base is zero and the 
exponent is not greater than zero, if the base is negative, or if 
arithmetic overflow occurs during computation. The result is 
calculated from 

base ** exponent = exp(exponent * log(base)). 

If the argument set is invalid, POS. INOEF. is returned and a 
diagnostic message is issued. If the argument set is valid, the 
computed result is returned to the calling program. 


3. ERROR ANALYSIS. 


The algorithm used in OTOX* is the same as that used in OTOX, • See 
section 3 of the description of routine DTOX. for an error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e' occurs in the base b and a small error e" 
occurs in the exponent p, the error in the result is given 
approximately by 

b**p * (p/b » e* ♦ log(b) * e"). 

The absolute error is approximately the absolute value of this 
expression. If the errors in the arguments are significant, the 
error in the result should be found by substitution of the possible 
argument values in the expression b p. 
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ROUTINE'S FUNCTION 


1.1. Type. fl FORTRAN exponentiating routine. It accepts an 
argument set comprising a double-precision quantity and a 
floating-point quantity* and returns a double-precision 
result« 

1.2. Purpose. To accent calls by value for OTOX. * calls generated 
by FORTRAN programs i*hich raise double-precision bases to 
floating-point exponents. 


?. METHOn. 

The input range is the collection of argument sets (b»p} for which: 
b is a definite in-range double-precision quantity, p is a definite 
in-range floating-point quantity, b is positive, if b is zero then p 
is greater than zero, and b^*p is in-range. 

The formula used is: 

b»»o = exp(o ♦ log b» 

where b > fl. Upon entry, OLNLOG. is called to compute log b, and o 
* log b is then computed in double-precision. OEXP* is called to 
compute expfp * log b), and the result is returned to the calling 
program • 


3. ERROR ANALYSIS. 

in,00n pairs (bfp) of random numbers were generated ( where b is 
doubIe-precIsi on and o is single-precision! with distribution the 
product of uniform distributions on (.5, 1.5) and CO, 1). The 
maximum absolute value of the relative error in the routine for 
these pairs was found to be 6.405 * 10**(-29l. 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e(b) occurs ir the base b and a small error e(o) 
occurs in the exponent p, the error in the result r is given 
approximately by 

r » (e(p) * log b p » e(b>/b). 
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1. ROUTINE'S FUNCTION. 


l.l. Type, ft FORTRAN exponentiating routine. It accepts an 
argueent set consisting of a double-precision argument and a 
complex argument, and returns a complex result. 

1.?. Purpose. To accept calls by name generated by FORTRAN 
programs raising double-precision quantities to complex 
quantities. 


2 , METHOO, 


If the base is real and the exponent is complex, then 
base ** exponent = x ♦ l.Y. 

where 

X = expfref exponent),Iog(base)).cos(iffl (exponent},logibase)) 

and 

y = exp(re(exponent).Iog(base)),s in(im(exponent)•Iog(base)). 
Upon entry the double-precision base is rounded to single-precision* 
and the resulting argument set is checked. The argument set is 
invalid if either number is infinite or indefinite, if the base is 
zero and the real part of the exponent is not positive, if the base 
is negative. if arithmetic overflow occurs during any stage of the 
computation, or if precision is lost through the arguments* being 
too large* If the argument set is invalid, a diagnostic message is 
issued and POS.INOEF, is returned. Otherwise, the result of the 
computation is returned to the calling program. 


3. ERROR ftNftLYSIS. 


The algorithm used in OTOZ* is the same as that used in OTOZ. . See 
the description of OTOZ. for an error analysis. 


<». EFFECT OF ARGUMENT ERROR. 


If e* and e** are small errors in the base b and exponent z 

respectively, then the corresponding error in b z is 

approximately ((z/b) » e* ♦ e** ♦ log (b)) » b**z. The absolute 

error will be approximately the absolute value of this. If e' or 
e" becomes significant, the error in the result should be 

calculated by substitution of the possible values of the arguments 
in the expression b z. 


100 


60498200 C 






1. ROUTINF'S FUNCTION. 

l.l. Type. A FORTRAN exponentiatirg routine. It accepts an 
argument set comprising a double-precision and a complex 
argument, and returns a complex result. 

1.?. Purpose. To accept calls by value for OTOZ# . calls generated 
by FORTRAN programs which raise double-precision bases to 
complex exponents. 


?. METHOD. 

The input range is the collection of argument sets (b*z) wherei b is 
a definite in-range double-precision quantity, z is a definite in¬ 
range comp 1exouantity. b is greater than zero, and b ** z and Jb ** 
z| are in-range. 

The formula used is* 

b**(u4-i*vl = explu^log b» * coslv*fog bJ 

«■ i . explu.log b) . siniv.log b) 

where b > 0 . Upon entry, the lower half of the double-precision 
base b is discarded, and ALOG. is called to compute log b . EXP. is 
called to compute exo (u.tog b). and COS=SIN is called to compute 
cosCv.log b) and sin (v.log bl. where u ♦ i.v is the exponent. The 
result is computed from the formula, and is returned to the calling 
orogram. 


?. ERROR ANALYSIS. 

10.nO’? pairs (b.z) (where b is double-precision and z is comolexi 
were generated with distribution the product of uniform 
distributions over (.5. 1.51 and 1-10.101 and (-2.pi.2.pi). The 
maximum modulus of the relative error in the routine was found to be 
5.605 * 10»»(-14) . 


4. EFFECT OF ARGUMENT ERROR. 


If a small error eCb) occurs in the base b and a 
occurs in the exponent z. the error in the 
approximately by 

w • (e(z) . log b z . e(b)^b> • 


smal I 
resul 


error 
t w is 


e (z) 
given 
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ROUTINE t ERF. 

1. ROUTINE'S FUNCTION, 

i.l. Tyoe, A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result, 

1.?, Purpose, To accept calls by value for ERF and ERFC from 
FORTRAN programs, ERF, computes the error function* 
ERFC, computes the comolementary error function, 1-ERF, 


METHOD, 


The input range is the collection of all definite floating-point 
quantities Cincludirg out-of-range values INF) except the range 
f P‘>,9227751502785<*» ♦■INF) for ERFC, which underflows. 


The routine calculates the smaller of erffabs(x)), erfcfabsix)), and 
uses the identities 

erff-x)s-erf(xl 
er ffxMl-erfclxl 


to compute the final value, which is the sum of a signed function 
and a constant. 

The forms used arei fysapslx)) 


£££ 


££££ 


\ 


I-INF,-5,625] 
(-S,F25,-,477) 
t-.477,0) 
to,*.4771 
(,477,5,625) 
15.625,8.0) 
t8.0,25,gi 
■♦INF 


-l.O 

-l,0ep2(y) 
-pi(y) 
♦pi (y) 
♦1.0-p2(y) 
♦ 1.0 
♦ 1.0 


♦ 1,0 


♦ 2.0 

♦2.0-p2f y) 
♦ 1. O^pKy) 
♦1.0-pl(y) 
P2(V) 
P2(y) 
under flow 
♦ 0.0 


where the constants .477 and 25,9 are inverse erffO.5) and inverse 
er f c ( 2 -®'^s) , which are approximately 0.47693627620447 and 
25,92277515027854. 


The function pi is a (5th order odd)/(8th order even) rational form. 
The functions d2,p3 are expf-x^)* (rational form), where d2 is (7th 
order)/(8th order) and p3 is (4th order)/(5th order). Since exp(- 
x*) is ill-conditioned for large x, expl-x^) is calculated by 
exo(uPS)=exDCu)+eps^exp(u ), where u=-x2 upper and eps=-x 2 lower. 

The coefficients for p2 and o3 are from Hart, Cheney, Lawson et al.. 
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7. EPROR ANALYSIS. 


The large error in d 2 and p3 is due to the large size of the 
rational forms and the additional error in expC-x®!. The 
polynomials in o2 and p3* while stable^ do not enjoy the high 
accuracy of most exponential-type approximations^ »«hlch» when 
evaluated using Horner's rule» sum the smallest terms first. 
Inverting x and reversing coefficients does not help due to the high 
error in divide. 

The maximum error in the approximations pi* p2# p3* scaled by 10*^®* 
is: 



&1 

Dl 

Xii 

rational form 

1.1 

4.9 

1.7 

coefficient rounding 

0.5 

0.8 

1.4 

round-o ff 

14.7 

110 

68 

upper bound 

16.3 

116 

71 

maximum observed 

12. 8 

27.9 

28.3 


In regions where a constant is added, that constant dominates and 
the error Is less than that shown. 


*». EFFECT OF ARGUMENT ERROR. 


For small errors in the argument x, the amplification of absolute 
error is (2/sqrtfpi))*expf-x*) and that of relative error is <2/ 
sgrt(piH*x*exp<-x*)/f <x> where f is erf or erfc. The relative 
error is attenuated for ERF everywhere and for ERFC of x<0.53. For 
x>n.53 the relative error for ERFC is amplified by approximately 2x. 

If the value of x is Known to more than single precision, the 
following sequence of FORTRAN may be used to compute a good value of 
ERFC when x is large! 

DOUBLE X 

DATA SQRTPI /<2/sqrtfpil>/ 


(compute X) 

SNGLX^SNGLfXI 

SHSNGLX=SNGL(X-SNGL(Xn 

YsERFC(SNGLX)♦SHSNGLX+SQRTPI*EXP(-SNGLX**21 


(Y is ERFC(XII 
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1. ROUTINE'S FUNCTION, 


1.1. Type. A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.?. Purpose. To accept calls by name and by value for EXP from 
FORTRAN programs. EXP computes the exponential function. 


?. MFTHOn, 

The input range (o this routine is the collection of all definite 
in-range floating-point quantities lying in the interval (- 
675.84»741.67), Upon entry, the argument x is multiplied by 16.7 
log(e)?. in double-precision, and the integral (n) and fractional 
(u) parts computed. The range reduction formula used here is 
expCx) = ? »» Ix/log0) 

= (2 *» 1/16) fl6 ♦ X / log 2) 

= (2 1/16) *♦ n * (2 1/16) u. 

If n = 16 ♦ q *■ r where q and r are integers such that 0Sr<16, 

exp(x) is finally given by 

exp(x) = 2 »* q * (2 ** 1/16) r » <2 1/16) ♦♦ u . 

q will be added to the exponent of the result. (2 ** 1/16) ♦♦ r is 
obtained from a IooK-uo table, and 12 ** 1/16) u is obtained from 
the following approximation 

<2 ** 1/16) »» u = 

u ♦ 2 » u * CpfOO) ♦ plOl) * u2) 


folOO) ♦ u2| - u ♦ Ip(OO) ♦ p(01) * u®) 

where the constants are given by 

qCnO) = 20.8137711965230361973 * 256 
plOO) = 7.2135034108448192083 * 16 
D(01) = .057761135831801928 / 16 , 

This aporoximation is described in Hart, Cheney, Lawson et al., 
$Computer Approximations* (New York) 1968, John Wiley S, Sons, 
op. 96-104. 


3, ERROR ANALYSIS. 


The maximum absolute value of the error of approximation of the 
algorithm is 5.000 • IQ -17 over the interval C-(log2)/16, 
floq2)/16). A graph of the error of approximation in the algorithm 
is given in figure 9. An upper bound for the absolute value of the 
error due to machine round-off is 1.868 * 10 ** -14 over the 
interval t(-log 21/16,(fog 2)/16], Hence an upper bound on the 
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absolute value of the error In the routine over this interval is 
1,873 * 10 ** -14. A bound on the routine’s error for any given 
argument x may be obtained by employing the multiplication formula 
for exp 

expfx + y) = expfx) • expfyl • 

The maximum absolute value of the relative error of approximation of 
the algorithm over (-log2/16 • log2/16) is 4.838 • 10**-17 . An 
upper bound on the absolute value of the relative error due to 
machine round-off and truncation is 6.890 * 10 -15 over t(-log 
2)/16,(log 21/161. So an upper bound on the absolute value of the 
relative error is 6.938 * 10 ** -15 over the interval t(-log 21/ 
16, Clog 21/161. 

For 10000 arguments chosen randomly from the foifoxing intervals, 
the following statistics on relative error were observed. 


Interval 
from to 

Mean 

Standard 

Deviation 

Minimum 

Maximum 

-673. 

-1. 

741. 

1. 

-3.Q12E-16 

-3.1Q0E-16 

2.181E-15 

2.223E-15 

-6.887E-15 
-6.769E-15 

5.193E-15 
5.028E-15 


4. EFFFCT OF ARGUMENT ERROR. 


If a small error e" occurs in the argument, the error in the result 
y is given approximately by y * e* • 
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I ROUTINE t HYP. C SINH l COSH ■> 
1. ROUTINE'S FUNCTION, 


1.1. Type. FORTRAN externa I functions. The routine accepts a 
floating-point argument* and returns a floating-point result. 

1.’. Purpose, To accept calls from FTN compiled code for 
comoutation of the hyperbolic cosine and sine functions. 


?. MFTHOO. 


The input range is the collection of all definite in-range, 
floating-point quantities Iyirtg within the interval 

t-1071*!ogf ?J ,1071»1ogC2) 1 = [-742.36063037 97,742.36 06 30 3797]. 

The formulae used to compute sinhixl and cos(x) aret 

X = n*log 2 ♦ a, lal<l/2*log 2, 

coshlx) = 2»*fn-l» * cosh(a) t 2*»(-n-l) * cosh(a) t 2**(n-l) » 

sinh(a) - 2*»(-n-l) * sinhial 

sinh(x) = 2*»(n-l> ♦ sinh(a) ♦ 2»»t-n-ll * sinh(a) + 2**(n-l) * 

coshfal - 2**(-n-ll * cosh(a) 
coshlal = 1 ♦ d(a) 

sinh(a) = a*a^* (s (3)♦•aa»(s C5) tj^/f^-a®) ) J 

dfa) = a2*(l/2 + a®»(c(4l+a2*CcC6)ta2»(c(8)+a®)*c(10)n) 

where 

Sl3) = .16666666666693558 
s(5» = -.005972995665652368 
h = 1.031539921161 
a = 72.10374670722 
c(4» = .041666666666488081 
cl6> = .0013888888952318045 
c(8) = 89.75473897315022 
C(IO) = 2.763250805803 » 10*»-7 

In the following description of the algorithm used, (XI) = x = 
argument on entry* entry is at SINH, or COSH, ? and on exit, (X6) = 
resuIt. 

a. If lxl>1071»log(2), go to step J. 

b. u *• I XI 

V «- +0 if x>Q 
- 0 I f x< 0 

d, n [u/log?t,5] = nearest integer to u/log2 

w •- u - n’^log 2, where the right-hand expression is evaluated in 
double-precision. 
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e« s *• (s (3)+*<2Is 151+Ji/Ca-'W*I)) 

d w*(l/2+«a(cC4>+wa(ct6»+w*(cC8)+Ha)*cC10nJ) 
a •- llfd-s>*2**l-n-l» 
b d*s 

f. If COSH, entry, go to step h, 

o. c «- (1/4 ♦ (l/4+b))*2**ln-H ♦ (2**(n-3) *■ l2**(n-3) - a)) 
X6 •- c with the sign stored in v. 

Go to step i* 

h. c ll^b»*2**fn-l) * a 

XG c 

i» Return* 

I* If infinite or indefinite arguaent, go to step (• 

H* formalize argument, 
u •• 1 xt 
w efl if k 5 0 
•0 if x<(J 

If Ixl < 107l»log 2, go to step d. 

I, Initiate error processing. 


X5 +IN0. 

if 

X 

is 

indefinite. 





♦INF. 

if 

X 

is 

infinite or 

too 

big. 

and 

positive, or COSH . 

-INF. 

if 

X 

is 

infinite or 

too 

big. 

and 

negative, and SINH . 


n. Go to step i, 

3, ERROR ANALYSIS. 


The maximum absolute value of relative error in the approximation of 
sinh over r-log2/2»log2/21 is 1.282 * 10»*-15 and of cosh over 
I-log2/2,log2/2J is 2.421 ♦ 10*»-16 . Computed upper bounds on the 
absolute value of relative error due to machine error in the 
computation of sinh is 2.392 * 10**-14, and of cosh is 1.024 * 10»»- 
14. Hence, upper bounds on the absolute value of relative error in 
the routine is 2.520 * 10**-14 for sinh, and 1.048 * 10**-14 for 
cosh. Graphs of the relative errors in the alogrithms used to 
approximate sinh and cosh over t-log 2/2, log 2/21 are given in 
figures 17 and 18. 


4. EFFECT OF ARGUMENT ERROR. 


Graphs of 
sinh and 
18. Tf a 
error in 


the relative errors in the algorithms used to approximate 
cosh over (-log 2/2,log 2/21 are given in figures 17 and 
small error u occurs in the argument x , the resulting 
sinhix) is given approximately by cosh(xl»u , and the 
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resulting error in coshtx) is given aoproximately by sinhixl*u« If 
the error u is not seal I• the addition formulae for sinh and cosh 
should be used to find the resulting errort 

sinhfx^u) = sinhfxIcoshful+coshfxIsinhCuJ 
cosh(x+u) = coshCxIcoshful♦sinhfxlsinh(u) 


no 
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1. ROUTINE'S FUNCTION. 

1.1. Type. Auxiliary functions from the FORTRAN CointBon Library. 
The routine accepts a floating-point argument and returns two 
floating-point results. 

1.?. Purpose. To accept calls by value from CCOS* « CCOS. t CSTN* 
and CSIN. for incidental computation of cosh and sinh. 


P, MFTHOO. 

The input range is the collection of all definite in-range floating¬ 
point quantities which lie in the interval (-741.67f 741.67). 

The hyperbolic cosine is computed by 

cosh(x) = .5 . (exp(x) ♦ exp(-x)). 

If !xl > .2?* the hyperbolic sinh is computed by 
sinhCx) = .5 . (exp(x) - exp(-x)). 

For Ixl < .22» the RacLaurin series for sinh is truncated after the 
term x*/9! and the resulting polynomial is taken as approximationt 
SinhCx) = X ♦ x3/3! ♦ xS/5! ♦ x'f/7i *■ x»/9! 


3. ERROR ANALYSIS. 


The maximum absolute value of the error of approximation for cosh(x) 
is F.flOO. 10*»(-17) and for sirh(x) is 1.464 . 10»*(-15), over the 
interval (-log2» log2). See the description of EXP. for details 
concerning the error of approximation to exp. An upper bound for 
the error due to machine round-off and truncation is computation of 
the MacLaurin polynomial is 8.198 . 1C**(-16)* A graph of the error 
of approximation in the polynomial for sinh is given in figure 10. 
An upper bound for the routine's error in the computation of cosh(x) 
is 7,184 . 10**(-14) and in the computation of sinhlx) is 7.148 
. 10**(-14) over (-log2. log2). 


4, EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in the argument x* the resulting error in 
coshCx)is given approximately by sinh(x).e'. and the resulting error 
in SinhCx) is given approximately by cosh(x).e'. 
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EmiiiiiE-i-iiimi- 


i, POUTTNE'S FUNCTION. 


i.l. Type. A FORTRAN exponentiating routine. It accepts an 
argument set consisting of a fixed-point and a double¬ 
precision argument, and returns a double-precision result. 

1.^. Purpose. To accept calls by name for ITOD* generated by 
FORTRAN programs raising fixed-point bases to double-precision 
quant ities. 


?. METHOD. 


The computation uses 

base ** exponent = expfexponent * log(base)). 

Upon entry, the fixed-point argument is converted to double- 
precision and the resulting argument set is checked. The argument 
set is invalid if the base is zero and the exponent is not greater 
than zero, if the base is negative, if either argument is infinite 
or indefinite, or if floating overflow occurs during the 
computation. If the base is zero and the exponent is negative, NEG. 
INF, is returned. If the argument set is otherwise invalid, POS. 
INOFF, is returned. In a 11 cases, if the argument set is invalid, 
a diagnostic message is issued. If the argument set is valid, the 
result is computed and returned to the calling program. 


?f, ERROR ANALYSIS. 


The algorithm used in TTOO* is the same as that used in ITOD* . See 
the description of routine TTOO. for the error analysis. 


4. EFFECT OF ARGUMENT ERROR, 


If a small error occurs in the double precision exponent, the 
resulting error in the result is given approximately by multiplying 
the argument error by the result, and then by the natural logarithm 
of the base. Thus, if the result is large, the effect of an 
argument error will be large. If the error in the argument becomes 
significant, the error in the result should not be calculated by 
this rule, but should be calculated from the function values. 
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i. ROUTINE'S FUNCTION, 


1*1« Tyoe. A FORTRAN exponentiating routine* It accepts an 
argument set comprising a fixed-point argument and a double¬ 
precision argument, and returns a double-precision result* 

1*?. Purpose* To accent calls by value for ITOO* , calls generated 
by FORTRAN programs which raise fixed-point bases to double- 
precision exponents* 


?. METHOn. 

The input range is the collection of all argument sets Cb»p) where* 
b Xs a definite in-range fixed-point quantity, p is a definite in¬ 
range double-precision quantity, bis greater than zero, and b**p is 
in-range* Upon entry b is floated, normalized and converted to 
doubie-precision* 

The formula used to compute the result is 
b**p = expfp • Iog b) • 

OLOG* is called to compute log b, then p.log b is computed in 
double-precision . OEXP* is called to compute expip.log b), and the 
result is returned to the calling program. 


ERROR ANALYSIS* 

10,0013 random argument sets fb,p) were generated, with distribution 
the product of a discrete uniform distribution over the integers 
1,1’,•*.,9 and a uniform distribution over (-1,1). The relative 
error in the routine was computed for each of the argument sets* 
The maximum absolute value of the relative error in the routine was 
found to be 2.466 * 10»*C-28). 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e* occurs in the exponent, the error in the result 
r is given approximately by r*e**1og b, where b is the base. 
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1. ROUTINE'S FUNCTION, 


1,1. Type, A FORTRAN exponentiating routine. It accepts an 
argument set consisting of t**o fixed-point arguments, and 
returns a fixed-point result. 


1 , 2 , 


Purpose, To accept calls by name for ITOJ» from FORTRAN 
programs which raise fixed-point quantities to fixed-point 
exponents. 


2. METHOD, 


Let b be the base and p the exponent. If p has binary 
representation 000,,,,Qfl0 i (n) i fn-1), , * i (i)i(0) where each ifj) 
(0<|<nl is 0 or 1, then 

D = + ifl),2i ifni,? n, 

and 

n = tlog(2)ol = greatest integer not exceeding log{2)p. 

Then 

b *» = Prod Cb»*C2»*n » 0<|<n i Kp = 1>. 

The numbers 1 - b®, b= b^, b«, b^,.,,, b»* (2 ** Clog(2)p]) are 

generated during the computation by successive squarings, and the 
coefficients i fH),, •,, i fn) are generated by sign tests of successive 
right shifts of p within the computer, A running product is formed 
during the computation, so that smaller powers of b may be 

discarded. The computation then becomes an iteration of the 

algor ithm i 

b o—bifp-l 

=(b,b) ( 0 / 2 ) if p is even 

= (b,b) ** ((o-l)/2),b if p is odd. 

Upon entry, the base is converted to floating-point, and the result 
of the computation will be later converted to fixed-point for 
return. The argument set is invalid if the base is zero and the 
exponent Is zero or negative, or if Integer overflow occurs during 
the computation. If the argument set is invalid, zero is returned 
and a diagnostic message is issued. If the base is non-zero and the 
exponent is negative, 1, -1 or 0 will be returned according as the 
base is 1 (or -1 with the exponent even), -1 (with the exponent 
Oddi, or other. The result of the computation is returned to the 
calling program. 


7, ERROR ANALYSIS - not applicable. 

4. EFFECT CF ARGUMENT ERROR - not applicable. 
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1. l?OUTTNE*S FUNCTION. 


l.l. Type. ft FORTRftN e«ponentiating routine. It accepts two 
fixed-point arguwents, and returns a fixed-point result. 

1.2* Purpose. To accept calls by value for TTOJ. generated by 
FORTRAN programs which raise fixed-point quantities to fixed- 
point exponents. 


2. METHOD. 


Special case 


Q 


Q 

s 

error 





J 

s 

error 

i f 

J<0 

-0 

** 

1 

s 

*0 



1 

** 

J 

= 

1 



-1 

e* 

J 

s 

♦ 1 or 

-IIJ even or odd) 

T 

** 

0 

s 

1 



T 

** 

J 

s 

0 if 

J<0 


I 

** 

2 

s 

I*I 



I 

** 

J 

s 

error 

i f 

T>2 and J>64 

T 

♦ ♦ 

J 

= 

error 

if 

I>2i« and J>3 


Let b be the base and p <>01 the exponent. If p has binary 
representation QOO•••OOi(n)i<n-l)*•*i (1)1(0) where each i(J)(0<J<n) 
is 0 or 1, then 

p = i<0).2o ♦ id) .21 ♦ i<2).22 i<n).2 n 

While p is even do 

b = b*f p = p/’2. 

Let r = b. 

While p > 1 do 
r = r*, 

if p is odd then r * r ♦ b. 
p = p/2. 

Now r contains the result. Floating point was used for r so that 
the remaining overflows could be caught by looking at the final 
exponent. 


ERROR ANALYSIS - not applicable. 


* 1 . EFFECT OF ARGUMENT ERROR - not applicable. 
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INE tITQX* 
ROUTINE’S FUNCTION 


1.1. Type. A FORTRAN exponentiating routine. It accepts an 
argument set consisting of a fixed-point base and a floating¬ 
point exponent, and returns a floating-point result. 

1.3. Purpose. To accept calls by name for ITOX* from FORTRAN 
programs which are generatec when fixed-point bases are raised 
to floating-point exponents. 


METHOD, 


Upon entry, the base is converted to floating-point, and the 
argument set is checked. The argument set is invalid if either 
argument is infinite or indefinite, if the base is negative, if the 
base is ^ero and the exponent is not greater than zero, or if 
floating overflow occurs during the calculation. If the base is 
zero and the exponent is negative, or if floating overflow occurs, 
POS, INF. is returned. If the argument set is otherwise invalid, 
POS.INDEF. is returned. In any case, if the argument set is 
invalid, an appropriate diagnostic message is issued. If the 
argument set is valid, the result is returned to the calling 
program. 


ERROR ANALYSIS, 


The algorithm used in ITOX» is the same as that used in ITOX, . See 
the description of ITOX, for an error analysis. 


EFFECT OF ARGUMENT ERROR. 


If a small error occurs in the floating-point exponent, the error in 
the result is giver approximately by multiplying the argument error 
by the result and then by the natural logarithm of the base. Thus 
if the result is large, the effect of an error in the exponent will 
be large. 
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ROUTINE'S FUNCTION 


1.1. Type. A FORTRAN exponentiation routine. It accepts an 
argument set fn,xi comprising a fixed-point argument n and a 
floating-point argument x, and returns a floating-point 
resuIt. 

1.?. Purpose. To accept caffs by value for ITOX. « calls which are 
generated by FORTRAN programs which raise fixed-point bases to 
floating-point exponents. 


?. METHOn, 

The input range is the collection of all argument sets tn.xf such 
that n is afixed-point quantity, x is a definite in-range floating¬ 
point quantity, x is positive and non-zero whenever n is zero, and 
n**x is in-range. 

The formula used is* 

n**x “ exp(x ♦ log n) , 
where n > 1. 

Upon entry, n is packed and normalized. Zero is returned if the 
base is zero. Otherwise, ALOG. is called to compute log n, and EXP# 
is called to compute exp ( x . log n f . The result is returned to 
the calling program. 


-T. ERROR ANALYSIS. 

50*[5,n00 pairs fn,xj of random numbers were generated with 
distribution the product of a discrete form of the right half of a 
Cauchy distribution, and a Cauchy distribution, n**x was computed 
for each of these pairs, first usingthe routine, and then using the 
double-precision routine. The maximum absoluteval ue of the relative 
error in the routine was 3.S29 * 10**(-12) for the 500,000 pairs. 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e* occurs in the exponent x, the error in the 
result r is given approximately by r * e* * log n, where n is the 
base. 
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gOUTlNE_»■ITQZ*■ 

1. ROUTINE'S FUNCTION 


1.1. Type. fl FORTRAN exponent!atlrg routine. It accepts an 
argument set comprising a real and a compfex argument, and 
returns a floating-point result. 

1.?. Purpose. To accept calfs by name for ITOZ* generated by 
FORTRAN programs which raise fixed-point bases to compfex 
exponents. 


?. METHOD. 


If n is a positive integer, and x and y are real, then 
n (x ♦ i.y) = expfx.logfnl) .cosfy. logfnH 

*■ i . exp (x. I ogfn I). s in (y. 1 ogCn) f 

Upon entry, the argument set is checked. It is invalid if the first 
argument is negative, or zero, if either argument is infinite or 
iniefinite, or if floating overflow occurs during the calculation, 
or if x^'log r is greater than 7A1.67. If the argument set is 
invalid, then a diagnostic message is issued, and POS. INDEF, is 
returned. Otherwise, the computation proceeds as outlined above and 
the result is returned to the calling program. 


T, ERROR ANALYSIS. 


The algorithm used in ITOZ* is the same as that used in ITOZ. • See 
the description of ITOZ. for an error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


Tf a small error occurs in the argument, the error in the result is 
given approximate?v by the product of the argument error, the result 
and the natural logarithm of the base. The absolute value of the 
error in the lesult will be given approximately by the product of 
the corresponding absolute values. If the argument error is 
significant, the error in the result should be found from 
substitution of the oossible argument values in the function. 
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I 


1. ROUTTNE-S FUNCTION. 


1.1. Type. A FORTRAN exponentiation routine. It accepts an 
argument set comprising a fixed-point quantity and a complex 
quantity, and returns a complex result. 

1.2. Purpose. To accept calls by value tor ITOZ. , generated by 
FORTRAN programs which raise fixed-point bases to complex 
exponents. 


?. MFTHOO. 

The input range to this routine is the collection of all argument 
sets (n,zl comprising a fixed-point quantity n and a complex 
quantity z such that z is definite and in-range, and such that* if n 
is zero then z is a positive non-zero real, im(z) . log n does not 
exceed pi.?^® (where n>0 and imfzl is the imaginary part of z), and 
the real number n re(z) is in-range. 

Upon entry, the fixed-point argument is packed and normalized, and 
then routine XTOZ* is called at entry XTOZ. to compute the result. 
The result is returned to the calling program. 


». ERROR ANALYSIS. 


100,000 pairs (n,z) of random numbers were generated with 
distribution the product of a discrete form of the right half of a 
Cauchy distribution, and the product of two Cauchy distributions. 
n*®z was computed for each of these pairs, first using the routine, 
and second using double-precision operations. The maximum absolut 
value of the relative error in the routine was found to be 3.054 
10**<-101 for these pairs. 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e(z) = e(x) ® i.e(y> occurs in the exponent z, the 
error in the result w is given approximately by w * log n * e(zl . 
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1 . 


ROUTINE'S function. 


1.1. Type. A FORTRAN external function. It accepts a dumffly 
argument and returns a floating-point result. 

1.?. Purpose. To accept calls by name for RANF and/ rANGET from 
FORTRAN programs. RANF computes pseudo-random numbers. 


?. NETHOO. 


RANF uses the multiplicative congruent,iai method modulo 2^®» i.e. 
x(nfl» = a * x(n) (mod 2^®) 

The library holds a random seed RANDOM, and a multiplier RANMLT. , 
The random seed can be changed to any value prior to calling RANF by 
use of the routine RANSET • Upon entry at RANF « the random seed is 
multiplied by the multiplier to generate a 96 bit product, and the 
lower A8 bits become the new random seed and is used to generate 
subsequent random numbers. RANDOM, has a default initial value of 
1T17 1274 T ?14 7741 31553 <241^®3 (nod • This new random seed 
is normalized and is returned as the random number. 

The multiplier RANMLT. is constant, and has a value of 2000 1207 
2642 7173 05658 . This multiplier can be shown to pass the Coveyou- 
Macpherson test as well as other statistical tests tor randomness, 
including the auto-correlation test with lag<100 and the pair 

triplet test (Referencet 0. E. Knuth, lijs_AcJi_gX_ Compu ter 

Programming , vol . 21. 

If RANF is called by name at entry point RANGET . the current seed 
of the random number generator is returned in the variable whose 
address is in XI . 


3. ERROR analysis - not applicable. 

4. EFFECT OF ARGUMENT ERROR - not applicable. 
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ROUTINE t ramse:i_ 

1. ROUTINE'S FUNCTION. 

1.1. Type. A FORTRAN external routine. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.?. Purpose. To accept calls by name for RANSET from FORTRAN 
programs. RANSET resets the seed of the random number 
generator. 


?. METHOD. 


The call supplied the new address of a Csuggested) new seed value in 
XI . If the new seed is fl. » the new seed value is made 
17171?7**321*i7741315FB (= .170R98394C4402317200I . Otherwise, the 

coefficient of the new seed is made odd if necessary Cby adding IB), 
and the exponent of the new seed value is set equal to -48 
C1717(8)). 


3. ERROR ANALYSIS - not applicable. 


4, EFFECT OF ARGUMENT ERROR - not applicable. 
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ROUTINE t SIMCQS, 

1. ROUTINE'S FUNCTION. 

1.1. Typ«. A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.?. Purpose. To accept calls by name and by value for SIN tat 
entry points SIN and SIN. respective!ylf and to accept calls 
by name and by value for COS (at entry points COS and 
COS. respectively). SINCOS. computes the trigonometric sine 
and cosine functions. 


?. HETHOO. 


The input range to this routine is the collection of all definite 
in-ranoe normalized floating-point quantities whose absolute values 
do not exceed pi ♦ 

Uoon entry, the range reduction 
y = 2/pi»x - n 

is performed in double-precision* where x is the argument* and n is 
an integer* and y is in I-l/?* 1/?1. Depending upon the sign of x 
and ntmod 4)* the result will be complemented or not* and a 
polynomial approximation (p(y) or qty)) will be chosen to give the 
result. The polynomial approximations pty) and qCy) are 

oty) = pi/?»y - y3»(s(0) ♦ s(l)*y2 ♦ s(?)*y** ♦ s(3)*y« ♦ s(4)*y»)2 
and 

qtyl = 1 - ya»(cC0) ♦ c(l)»y2 ♦ ctEl^y** ♦ c(3)*y« ♦ c(4)*y*)2 

The coefficients are 

StO) = 8.037l8gi69T670e * 10*»-1 
stlJ = -4,95774235001375 * 10**-2 
St?) = 1.38346449783347 * 10**-3 
s(3) = -1.44725130681196 ♦ 10**-S 
s(4) = 1.54733311005155 * 10**-7 
ctfl) = 1. 110720734539535 
ctl) = 1.1419139843400? * 10**-1 
ct?) = -3.521949713998275 * 10**-3 
c(3) = 5.172606069276518 * 10*»-5 
c(4S = -4.413?82528387igi » 10**-7. 

The polynomial approxinetiors pty) and qty) are minlmax 
apps^rximations to their cor c .pcsiding functions over I-pi/4, oi/4 1, 
(The algorithm and consiiC-^ ate copyright 1970 by Krzysztof 
Frankowski* Computer Inform, tic-f and Control Science* University of 
Minresota* 55455, and ar^ eppfoyed under licence. Coding is by 
Larry Liddiard* University of Minnesota.) 
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ERROR ANALYSIS 




A graoh of the error of approximation in the algorithm for sinfx) 
over t-pi/4» pi/41 is given in figure 3* and for cosCxl over [-pi/4» 
pi/41 in figure 4- The maxiwuiB absolute value of the error of 
approximation in the algorithm for sin(x) over [-pi/4» pi/41 is 
5.670 * and for cos(x» is 2.972 * 10**-15. Upper bounds 
for the error due to machine error in the computation of sin(x) and 
cosCx) were established at 2.898 * 10*»-14 and 4.444 * 10**-14 
respectively. Hence upper bounds on the error in the routine are 
2.955 * 10**-14 and 4.741 ♦ 10*»-14 for sinCx) and cosfxl, 
respectiwely. 

The maximum absolute value of the relative error of approximation in 
the algorithm for sinfxl over t-pi/4» pi/41 is 4.098 * 10**-14 and 
for cosCxI is 6.285 * i(l**-14. Upper bounds for the absolute value 
of the relative error due to machine error in the computation of 
sinixl and cos(x) were established at 8.049 * 10**-16 and 4.204 * 
10**-15 respectively. Hence upper bounds on the absolute value of 
the relative error in the routine were established at 4.178 * io*»- 
14 and 6.70*=^ ♦ 10**-14 for sinixl and cos(x) respectively. 

For 1000 arguments chosen randomly from the following intervals for 
the entry points shownt the associated statistics on absolute or 
relative error were observed. 


Entry Error 
Point 


In tervaI 
from to 


Mean 


Standard Minimum Maximum 
Deviation 


COS. 


Re 1ative 

-.7854 

.7854 

-5933E-17 

1.596E-15 

-7.346E-15 

6.962E-15 

Abso lute 

-3.1416 

3.1416 

-7.524E-18 

1.317E-15 

-4.674E-15 

4.809E-15 


-1012 


8.138E-19 

1.248E-15 

-5.443E-15 

4.843E-1S 


SIN, 


Relative 

-.7854 

.7854 

3.035E-16 

1.984E-15 

-6.448E-15 

6.739E-15 

Absolete 

-3.1416 

3.1416 

-2.504E-18 

1.133E-15 

-5.648E-15 

5.174E-15 


-1012 

10^2 

-6.872E-18 

1.254E-15 

-4.187E-15 

5.353E-15 


4. EFFECT OF ARGUMENT ERROR. 


Tf a small error e* occurs in the argument x. the error in the 
result is given approximately by e* * cosfx) for sinixl* and -e* ♦ 
sinCx) for cosix) . 
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ROUTINF t stnc<;d. _ 

1. ROUTINE'S FUNCTION. 


1»1» Type. A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.?. Purpose, To accept calls by value for SIND and COSO, the 
trigonometric sine and cosine function with argument in 
degrees. 


?. METHOD. 


Argument range* (-Z***, 

Routine OEGCOM, is called to subtract the necessary multiple of 90 
from the argument to put the result in t-45, ♦AS) and multiply the 
reduced value by pi/180. The appropriate sign is copied to the 
value of the appropriate function (sine, cosine) as determined by 
these identities* 

slr(X+360<») = sin(X) 
sin(Xil80O) = -sin(X) 
sin(X+90O) = cos(X) 
sin(X-qoo) a -cos(X) 
oos(Xi360®) = cos(X) 
cos(X+180®) = -cos(X) 
cos(X+go®) = -sin(X) 
cos(X-90O) a sin(X) 


T. ERROR ANALYSIS. 


The reduction to I-4S, *45) is exact? the constant pi/180 has 
relative error 1.37E-15, and the multiply by this constant has 
relative error S.33E-15, for a total error of 6.7E-15. Since errors 
in the argument of SIN and COS contribute only (pi/4) of t^^eir value 
to the result, the error due to the reduction and conversion is at 
most S,?6E-1S. The total error in SIND and COSO is at most this 
value plus the maximum error in SINCOS. on [-pi/4, ♦pi4), namely 
7,31E-15, for a total of 1E.57E-15. The maximum observed error in 
lOOClOn points in the interval [0,360) was 9.96E-15 for SINO and 
9.95E-15 for COSO. 


4. EFFECT OF APGUMENT ERROR. 


Errors in the argument X are amplified by X/tan(X) for SINO and 
X*tan(X) for COSO. These functions have a maximum value of pi/4 in 
[-45°, *45®1 but have poles at even (SINO) or odd (COSO) multiples 
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When 


of 99®* and are large between multiples of 90® If X is large. 
X is known to double orecisiort the following code may be used* 

FUNCTION SlNOOfXI 
DOUBLE X 

NTNT(X)«X*SIGNC0.5,XI 

K*0 

GO TO 1 
ENTIfY COSOO 
Ksl 

1 NsNINTfSNGLfXl/901 
?*X-N*90 

IfCK.NE.HOOflAeSfBI*?)} GO TO ? 

YsSINOCZ) 

60 TO 3 

2 YsCOSOfZI 

7 TF(K*2-l.E0,M00fN,?llY s-y 
IFCMOO flABSINI,41.GE.2IYS-Y 
SIN00*Y 
RETURN 
END 
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S2UlIIg-l.-SaBT, 

!• ROUTINES FUNCTION 

1»1« Type* A FORTRAN external function* It accepts a real 
argument and returns a real result* 

1*?* Purpose* To accept calls by reference for SORT from FORTRAN 
programs* SORT computes the square root function* 


?* METHOD* 


The argument is loaded into XI and tt>e call is converted to a 
SORT* call. 


T* ERROR ANALYSIS - see SORT* 


4. EFFECT OF ARGUMENT ERROR - see SORT* 
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ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function* It accepts a real 

argument and returns a real result. 

1.3. Purpose. To accept calls by value for SORT. from FORTRAN 
programs. SORT, computes the square root function. 


?. METHOD. 

The argument range is the set of all positive or zero floating point 
numbers. The identity 

sqrtCy*2ep)=sqrt(y)*2eCp/2l 

is used to reduce the range to C0.S*1) with p having an integral 
value. An initial approximation is made using one of eight linear 
approximations to sqrt on this interval* giving at least 12 bits of 
accuracy. Two Heron's rule iterations are made to obtain 49 bits. 

If p is even* the normal Heron's rule is used* 

compute xO* an approximation to x*sqrtfy) 

xl=0.5*(x0+y/x0) 

x2=0.5*(xltv/xll 

If p is odd* scaling is done between steps so as not to affect the 
accuracy of the final result* 

compute xO 
xl=0.5»<x0*y/x0J 
xl'=xl*sqrtC2) 
x2=0.«5»ixl'e(2*y) /xl') 

which accomplishes the multiply by 2ef1'2)=sqrt(2). 

The scaling by 2etp^2] Ctui denotes truncation! is done by packing 
the appropriate exponent with the coefficient of (2*x2l. The square 
root of a number one ulp below an even power of 2 is explicitly 
forced to one ulp below the square root of that power of 2 to make 
packing work* e.g.* sqrt <4-eps) would be 1.0 but is forced to 2- 
eps • 

The sqrt C2> scaling is fudged slightly so that the error is 
centered after this scaling* picking up one bit at that point* 


3. ERROR ANALYSIS. 


The maximum error in the initial approximation is .000218. Since 
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the effect of a Heron's iteration is to square and halve the 
relative errort the algorithm error is 7,08E-17, 

Round-off error is insignificant until the last Heron's rule step» 
which has the form x+y/x* where the quantities being summed are 
almost equal* Since the error in Heron's rule is always positive* x 
is too large* so y/x is too small* i.e.* x>y/x* The error in the 
divide is in (-7. tE-l?5* Ql and in the rounded odd is in tO*-^t,55E- 
15)* so the total round-off error is less than 3.55E-15 in absolute 
value, (Error in divide is halved because x=y7x approx*) 

The upper bound on relative error is then 3*62E-15. The maximum 
observed relative error in 100000 randomly chosen point in the 
interval £0,5*2) was 3*5gE-15* 


U, EFFECT OF ARGUMENT ERROR. 


For small error in the argument 
error is 1 /(2’^sqrt (y)) and that 


y the amplification of 
of relative error is 0*5* 


absolute 
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ROUTINE’S FUNCTION 


1,1. Type. An auxiliary routine. 

Purpose. To provide a link between routines in the math 
library, and the system error processor. 

METHOD. 

Execution proceeds as follows. 

a. Enter SYS=AIO and additionally save registers X3 and X4. 

h. Read up entry point SYSAIO. and store it at entry point 
SYSIST. . 

c. Long Jump to MORGUE. « 

See the method description of SYS=1ST for further details. 

3. ERROR ANALYSIS. 

Not applicable. 

4. EFFECT OF ARGUMENT ERROR. 

Not applicable. 
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eQUIINE.8 S YS-iST 


1. ROUTINE'S FUNCTION. 


1.1. Type. An auxiliary routine. 

1.?. Purpose. To provide a link between routines in the math 
Iibrary and a system error processor. 


?. METHOD. 


Execution proceeds as follows at MORGUE, t 

a. Enter SYSIST and save registers XI , X? , X6 and AO , 85 ♦ 36 

and 87. 

b. Read the return Jump word used to enter the routine which called 
SYS=1ST or SYS®AIO . If this word has the format* 

♦ RJ <entry point> 

VFO 3071 

then go to f. below. 

c. Read the communication cell SYSAIO* • Insert in its lower 18 
bits the address of the trace word in routine SYS=1ST • Store 
the result in cell RJERR which will be executed at step e. 

dl. Test the argument in the register Indicated by the contents of 
82 . Set X2 to the first word address of an error message as 
follows * 

Condition Message 

Infinite ARGUMENT INFINITE 

Indefinite ARGUMENT INDEFINITE 

Other ARGUMENT <partlal message from address supplied 

in 82 > 

Set XI to the error number, and AH to the first word address of 
the parameter list for nort-standard error recovery. 

e. Execute word RJERR • This will link the routine to the system 
error processor. 

f. Restore registers XI , X2 , AO » B5 , 36 , 37 . Move the 
entering contents of X6 into register X5 . 

g. Set X6 and X7 to elNO. . 

h. Return to the calling program. 


3. ERROR ANALYSIS. 


Not applicable. 


4* EFFECT OF ARGUMENT ERROR. 


Not applicable. 
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ROUTINE'S FUNCTION, 


1.1, Type. fl FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result, 

1,?, Purpose, To accept calls by value for TAN computes the 
trigonometric tangent. 


?, NAMES 

• 



2.1, 

I dent name - TAN 



2.2. 

UPDATE deck name 

- TAN 


2.3. 

Entry point name 

- TAN 


3. CALLS 

• 



3,1. 

Source of calls, 
compiled under 

mentioning TAN in 

From FTN compiled code mentioning 
control card option T. 0, or 
an EXTERNAL statement. 

TAN and 
OPT*0 or 

3,2, 

Format of calls, 
lump to TAN. 

Call by reference. Entry is made 

by return 

3.3. 

Format of return. 

The result is returned in X6, 



4, CALLED ROUTINES. 

TAN, at entry point TAN. to compute the result, 

F. HETHOO. 

2. HETHOO, 

The argument is loaded into XI and the call is converted to a 
TAN. call. 

3. ERROR ANALYSIS - see TAN, 

4. EFFECT OF ARGUMENT ERROR - see TAN. 
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1, ROUTINE*S FUNCTION. 


1.1. Type, fl FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.?. Purpose. To accept calls by name for TAN. from FORTRAN 
programs. TAN. computes the trigonometric tangent functions. 


?. METHOD. 


The input range is the collection of all definite. in-range 
floating-point quantities in the interval . 

The identities 

i) tan(xl = tan(x«-k*pi/2> k even 
ii) tan(xl=-l.0/tanlx+pi/2> 

are used in the form 

iiil tan(x)=tan<Cpi/21»fx*2/oi*kl> k even 
iv) tan(x>=-l.n/tan< (pl/2l*<x*2/pl + l)l 

to reduce the evaluation to the interval t-0.5*+0.51 using an 
approximation for tan<(pi/2»y). The reduction is done by 
multiplying x by 2/pl and subtracting the nearest integer, rounding 
the result to single. 

The function tanf Ipi/2)*yI is approximated with a rational form. 
(7th order oddl/(6th order even), which has minimax relative error 
on the interval I-F.5,+0.5). The rational form is normalized to 
make the last numerator coefficient (l*eos) where eps is chosen to 
minimize rounding error in the leading coefficients. 

If identity (iv) is used. i.e.. if the integer subtracted is odd. 
the result is negated and inverted by dividing -Q/P instead of P/Q. 


?!. ERROR ANALYSIS. 


The range reduction, the final add in each part of 
form, the final multiply in P and the divide dominate 
Each of these operations contributes directly to the 
and each is accurate to about 1/2 ulp (unit in the last 
maximum relative errors are 


the rational 
the error, 
final error, 
place). The 
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$aBiount*10*^s$ 


fsource of error* 


range reduction 3.6 

rational form *02 

coefficient rounding <.08 

round-off 14.2 

upper bound 18.0 

maximum observed 14.5 


4. EFFECT OF ARGUMENT ERROR. 


For small errors in the argument xt the amplification of absolute 
error is sec*(xl and that of relative error is x/Csin(x)*cos(x)1. 
which is at least 2x and may be arbitrarily large near a multiple of 
oi/2. If X is known to more than double precision* the tangent 
addition formula may be used if x is less than 3F7 i 

DOUBLE X 

Icompute X) 

T=TAN(SNGL(X)J 

F*SNGL(X-SNGLIXI 1 


YsT*S»Cl*T**2l/C1-S*T1 

CSaTANfSI if X<3E7|. This approximation may give less than single 
precision when S»T is near 1.0* where it is more accurate than 
TANISNGLIXll but less accurate than SNGL(OTAN(X>>. 
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ROUTINE t TAND. 

1. ROUTINERS FUNCTION 


1«1. Type.A FORTRAN externa! function. It accepts a floating-point 
argument and returns a floating-point result. 

1.?. Purpose. To accent calls by value for TANO, the trigonometric 
tangent function with argument in degrees. 


?. METHOh. 


Argument ranget *2^®} except odd multiples of 90. 

Routine OEGCOM. is called to subtract the necessary multiple of 90 
from the argument to out the result in t-45, *45) and multiply the 
reduced value by pi/18Q. Routine TAN. is called to compute the 
tangent, and the result is negated and inverted if the multiple was 
odd, using these Identities* 

tan(X+180O) = tan(X) 

tanCXi90O) s -l/tan<X) 


■’. ERROR ANALYSIS. 


The reduction to t-45, ♦45) is exact? the constant pi/180 has 
relative error 1.37E-15, and the multiply by this constant has 
relative error 5.33E-15, for a total error of 6.7E-15. Since errors 
in the argument of TAN are amplified by at most oi/2, the error due 
to reduction and conversion is at most 10.52E-15. The error in the 
final divide is at most 7.11E-15, and the error in TAN. is at most 
14.54E-15, so an uoper bound on error in TAND is 32.17E-15. The 
maximum observed error in 100000 points in the interval 10,360) was 
17.72F-15. 


4. EFFECT OF ARGUMENT ERROR. 


Errors in the argument X are amplified by at most X/(sinfX)*cosCX)). 
This function has a maximum of pi/2 within £-45° ♦45®] but has poles 
at all multiples of 90® except 7ero and is at least 2»X elsewhere. 
When X is known to double precision and one of the above problems 
exists, the following code may be used* 

fcompute X in double) 

N=NINT fSNGL fX ) /95) 

Y=TAN0(SNGL(X-N*90)) 
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IF {N00{N,2) .EQ.QI GO TO 1 
IF(Y,FQ.0J <error> 

Y=-1.0/Y 
1 CONTINUE 

which always returns an accurate value since the 
exact. 

CNoteS NINTIX) = IFIX (X •^SIGNI0 .5, XJ » , the nearest 


range reduction Is 
Integer*) 
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i. ROUTINF'S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a floating 
-point argument and returns a floating-point result. 

1.?. Purpose. To accept calls by value for TANH from FORTRAN 
programs. TANH computes the hyperbolic tangent function. 


P, METHOD. 


The argument is loaded into XI and the call is converted to a 
TANH. call. 

T. ERROR ANALYSIS - see TANH. 

4, EFFECT OF ARGUMENT ERROR - see TANH, 
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I ROUTINE t TANH, 

1. ROUTINE'S FUNCTION. 


1.1. Type. A FORTRAN external function. It accepts a floating¬ 
point argument and returns a floating-point result. 

1.2. Purpose. To accept calls by value for computation of the 
hyperbolic tangent, including converted calls from TANH . 


2. METHOD. 


The input range is the collection of all definite floating-point 
quantities in the range I-INF,*INF1, 

The identity tanhi-x)=-tanh<xl is used to reduce the range to 
ti3,tINFl. For abs(x)>17.50, the best machine representation of 
tanhCxJ is sign(l.Oyx), so the range is further reduced to 
10,17.50) . 

The identities 

tanh(x)=p(xl/a(x) approximately, on £0,0.551 
tanh(x)=:l-2/(exp(2*xl <■!) 
exp(2*xl = (1 + t anh fX)I/(l-tanh fx)} 
expf2*x)=?'tn»exp {2*fx-n*l n<2)/2)) 

may be combined to get 

tanh(x) = l-2*(q-ol / (f a-p)>2en*<q>pn 

where n is chosen to be nintCx*2/In(21) and o,q are evaluated on x- 
n*lnl?)/2. This choice of n minimizes absCx-n*ln(2)/2I. 

When x<0.55 the approximation p<x)/q(x) is used. Since 
tanhlx<0.55)<0.5, the form 1-r would suffer from cancellation in 
this range. 

The approximation p/q is a minimax Crelative error) rational form 
i.e., (5th order odci)/(6th order even). The coefficients are scale 
so that (x*2/In(2)-n) may be used Instead of (x-n*ln(2)^2), 

simplifying the range reduction. The coefficients are further 
scaled by an amount sufficient to reduce truncation error in the 
leading coefficients without otherwise affecting accuracy. 


7. ERROR ANALYSIS. 


The algorithm error (due to finite approximation, coefficient 
truncation) is 1.7E-15. For abs(x)<0.55 the form o(x)/qfxl is used, 
and the final operations z«-x*2/ln(2) and tanh«-(z*(o0+sfflal 11) / 
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CaO«-small) dominate the error. The upper bound on the error here Is 
18.0E-15* the maxleum observed was 13.0E-15. 

For abs(x)>1.25 the final subtract* l.C-smallt dominates and an 
upper bound on the error is 4«2E-15t the fflaxiraum observed was 3.8E- 


For 0.55<absixl<1»25 the final operation is 1-R where R becomes 
smaller as x approaches 1.25, so the worst relative error is near 
8.55* ramely (contribution from Rl+Cerror in final sum)* where 
R=2*(q-p)/((q-p)*4*(q^pl). An upper bounds 16.7E-15* maximum 
observed: lO.OE-15. 

Relative Error: 

$source of error* Serror 

rational form 0.5 
coefficient rounding 1.2 
round-off 16.5 
upper bound 18.2 
maximum observed 13.0 


4. EFFECT OF ARGUMENT ERROR. 


For small errors in the argument x* the amplification of the 
absolute error is l/cosh2(xl and of relative error is x/ 

(sinh(xl*cosh(xl). Both have maximum values of 1.0 at 0 and 

approach 0 as x gets large. 
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EfiUTINE IJilQIll 


i. ROUTINE’S FUNCTION, 


1.1. Tyoe, fl FORTRAN exponentiating routine. It accepts an 
argument set comprising a floating-point and a double¬ 
precision argument, and returns a double-precision result, 

1.2. Purpose, To accept calls by name for XTOO* generated by 
FORTRAN programs whicH raise floating-point bases to double¬ 
precision exponents. 


2. METHOD. 


The formula used is* 

base ** exDonent = expCexponent » iog(basel). 

Upon entry, the argument set is checked. It is invalid if either 
argument is infinite or indefinite, if the base is negative, if the 
base is zero and the exponent is not greater than zero, or if 
floating overflow will occur during the computation. If the 
argument set is invalid, a diagnostic message is issued and 
POS.TNOFF, is returned. If the argument set is valid, the result 
is returned to the calling program. 


2, EP20R ANALYSIS. 


The algorithm used in XTOO* is the same as that used in XTOO, • See 
the description of routine XTOO, for ar error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


If a small error e* occurs in the base 
occurs in the exponent p, the error in 
aporoximate Iy by 

b**p * fp/b * e* ♦ logfb) * e**l , 

The absolute error is approximately the 
expression. If the errors in the argument 
error in the result should be found by subst 
argument values in the expression b *♦ p. 


b and a small error e*" 
the result is given 


absolute value of this 
are significant, the 
itution of ♦he possible 
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1. ROUTINE'S FUNCTION, 


1,1. Type, A FORTRAN external function. It accepts an argument 
set comprising a floating-point and a double-precision 
argument, and returns a double-precision result, 

1,?. Purpose, To accept calls by value for XTOO, , calls generated 
by FORTRAN programs which raise floating-point bases to 
double-precision exponents. 


?, METHOD, 

The input range is the collection of argument sets lb,p) where! b is 
a definite in-range floating-point quantity, p is a definite in¬ 
range double-precision quantity, b is greater than zero, and b**o is 
in-range. The result is computed according to b**p = exp(p,log b), 
where b is converted to double-precision upon entry, and all 
operations are carried out in double-precision. The result is 
returned to the calling program. 


?, ERROR ANALYSIS, 

10,00n argument sets Cb,pl were randomly generated, with 
distribution a product of uniform distribution on f,?,l,?} and (- 
10,101, The relative error in the routine was computed for each of 
the argument sets. The maximum absolute value of the relative error 
was found to be 1,163 • 10**C-25), 


%. EFFECT OF ARGUMENT ERROR. 


If a small error 
occurs in the 
approximately by 
r . lelpl 


elbJ occurs in the base b and a small error 
exponent p, the error in the result r is 


log b ♦ p , efbl/b) 


elpl 

given 
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ROUTINE » XTOI* 

1. ROUTINF-S FUNCTION. 


1.1. Type. ft FORTRAN exponentiating routine. It accepts an 
argument set consisting of a floating-point and a fixed-point 
argument, and returns a floating-point result. 

1.2. Purpose. To accept calls by name for XTOI* generated by 
FORTRftN programs which raise floating-point bases to fixed- 
point exponents. 


2. METHOO. 


Load arguments and call XTOI. 


?. ERROR ANALYSIS. 


Not applicable, since the only errors are round-off errors. See the 
description of XTOI. . 


4. effect of ARGUMENT ERROR. 


If a small error e* occurs in the base o. the error in the result is 
given approximately by 

b**(p«i) * p ♦ e*. where p is the exponent. 

If the error in the base becomes significant, the error in the 
result must be found from substitution of the possible values of the 
base b into the expression b •• p. 
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1, ROUTINE'S FUNCTION, 


1.1, Type, A FORTRAN exponentiating routine. It accepts an 
argument set consisting of a floating-point quantity and a 
fixed-point quantity, and returns a floating-point result, 

1.2, Purpose, To accept calls by value for XTOI# , generated by 
FORTRAN programs which raise floating-point quantities to 
fixed-point exponents. 


?, METHOD. 


Special case 

X indefinite t error 
X infinite ! error 
0, 0 * error 

X *♦ p = 1.0 

(number of bits in I)*-(number of bits in scale of Xl>8 * use 
careful code 

X *♦ I = 1,0 /(X **(-TH if KO 

Quick version? Walk through the binary representation of I, 
starting with the most significant bit. For each bit, square the 
result (which was initialized to Xl? if the next bit is on, also 
multiply by X. 

Careful version? Scale X to be between 0,75 and 1,5, remembering the 
exponent. Walk through 10 bits of I in the quick-vers ion way. Then 
repeat (sea te^-wa Ik) until I is used up. Invert if necessary. 
Carefully decide if the exponent is too big. If OK, exit. 


:t, ERRCR analysis - not applicable. 


4. EFFECT OF ARGUMENT ERROR. 

If a smalI error e" occurs in the base b, then the error in the 
result is given approximately by p * b**(p-ll » e* , where p is the 
exponent. If the error e* becomes significant, we can only say that 
the absolute error in the result is bounded above by 
Ipl ♦ maxdbl , lb *■ e*l))*»(p-l) ♦ le'l . 
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ROUTINE*'S FUNCTION 


1.1. Type. A FORTRAN exponentiating routine. It accepts an 
argument set consisting of two floating-point arguments and 
returns a floating-point result. 

1.?. Purpose. To accept calls by name for XTOY* generated by 
FORTRAN programs whicti raise floating-point bases to floating¬ 
point exponents. 


?. METHOD. 


The formula used is 

base “f-* exponent = exp (exponent . 1 ogCbase'* > . 

The argument set is checked upon entry. It is invalid if either 
base or expcnent is infinite or indefinite* if the base is negative, 
if the base is zero and the exponent is not greal^r than zero, or if 
floating overflow occurs during the computation. If the argument 
set is invalid. POS.INOEF. is returned and a diagnostic message is 
issued. Otherwise, the result of the computation is returned. 


3. ERROR ANALYSIS. 


The algorithm used in XTOY* is the same as that used In XTOY. • See 
the description of routine XTOY. for an ennon analysis. 


L. EFFECT OF ARGUMENT ERROR. 


If a small ennon e" occurs In the base b and a small error e*** 
occurs in the exponent p. the error in the result is given 
approximately by 

b»*o » (p/b * e* ♦ log(bl * e**). 

The absolute error is approximately the absolute value of this 
expression. If the errors in the arguments are significant, the 
error in the result should be found by substitution of the possible 
argument values in the expression b ** p. 
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1. ROUTINE'S FUNCTION, 


1,1, Type, fl FORTRAN exponentiation routine. It accepts an 
argument set comprising two floating-point arguments, and 
returns a floating-point result. 

1,?, Purpose, To accept calls by value for XTOY, , generated by 
FORTRAN programs which raise floating-point bases to floating¬ 
point exponents. 


P. METHOD. 

The input range is the collection of all argument sets <b,e) for 
whicht b and e are definite in-range floating-point quantities, b is 
positive and non-zero, and b»*e is in-range. 

The formula used is* 

b»*p = expIp * log b), 
where b > 0, 

Upon entry, ALOG, computes log b, and then EXP, computes 
expfp ♦ log b), 

The result is returned. 


T, ERROR ANALYSIS. 

*500,011(5 pairs (b,p) of random numbers were generated with 
distribution the product of the right half of a Cauchy distribution, 
and a Cauchy distribution,b ** p was computed for each of the pairs, 
first using the routine, and then using the double-precision 
routine. The maximum absolute value of the error in the routine was 
4.58;? ** 10»*(-12) for these 500,000 pairs. 


4. effect of ARGUMENT ERROR, 

If a small error e(bl occurs in the base b, and a small error e(p) 
occurs in the exponent p, the error in the result r is given 
approximately by 

r * (log b * e»*p ♦ p » (efbM/b) , 
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ROUTINE 


1, ROUTTNE-S FUNCTION. 


1.1. Tyoe. A FORTRAN exponentiating routine. It accepts an 
argument set comprising a floating-point base and a complex 
exponent, and returns a complex result. 

1.?. Purpose. To accept calls by name for XTOZ* generated by 
FORTRAN programs raising floating-point quantities to complex 
exponents. 


?, METHOD. 


If the base b is real and the exponent z = x ♦ i * y where x and y 
are real, then 

b ** z^u+i^v. 

where 

u = exp (x ♦ logfbll * cos (y ♦ log(bl) 

and 

V = exp fx * log(b)l * sin Cy ♦ logCbl). 

ALOG.. EXP. and COS=SIN are called to evaluate these expressions. 
The argument set is checked upon entry. It is invalid if either 
base or exponent is Infinite or indefinite. if the base b is 
negative, if the base is zero and the real part of exponent z is 
greater than zero, if y * log fbl is so large that precision is lost 
in the computation, or if floating overflow occurs during the 
computation. If the base b is zero, y is zero and x is less than 
zero. POS, INF. is returned. If the argument set is otherwise 
invalid. PCS, INOEF, is returned. In either case, a diagnostic 
message is issued. If the argument set is valid. ALOG. . EXP. and 
C0S=5iIN are called during computation. The result is returned to 
the calling program. 


T. ERROR ANALYSIS. 


The algorithm used in XTOZ* is the same as that used in XTOZ. . See 
the description of routine XTOZ. for an error analysis. 


4. EFFECT OF ARGUMENT ERROR. 


If a small error eCbl occurs 
and e(yl occur in the real 
frespectively) of the exponent 
is given approximately by 
e(r) = b»»z ♦ logfb»*z*(<e <x» ♦ 
The absolute error in the result 


in the base b. and small errors e(x) 
and imaginary parts x and y 
z. then the error e(r) in the result 

i*efyn/z ♦ e (b)/<b*log(b))) . 
is approximately the absolute value 
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of this expression. If the error in an argument becomes 
significant, the error in the result should be found from 
substitution of possible argument values in the expression b ** z. 


152 


60498200 C 



SQUTIN^.I.xjQZx 

1. ROUTINE'S FUNCTION 


1.1. Type. A FORTRAN exponentiating routine. It accepts an 
argument set CxtZl where x is a floating-point quantity and z 
a complex quantity, and returns a complex result. 

1.2. Purpose. To accept calls by value for XTOZ., calls which are 
generated by FORTRAN programs which raise floating-point bases 
to complex exponents. 


2. METHOD. 

The input range is the collection of all argument sets 
<x,zl ( = X, u *■ i*vl) 

such thati X is positive, if x is zero then u = 0 and v is positive 
and non-zero, both x and z are definite and in-range, floating over¬ 
flow does not occur during the computation of x ** u fi.e., 
{u.log(x)) < 741.67 , and Iv.togfxH < pi.2^® . 

The formula used ist 

x*’^(u+i*w) = e**(u*log<xl * cos{v®log(x) I 

* i ♦ e**(u» logfx) * sinCv*I og(x) >. 

Upon entry, the base is checked. If it is zero, zero is immediately 
returned to the calling program. Otherwise, ALOG. is called for 
computation of log x, and then COS=SIN is called for computation of 
coslv.logfx)) and sin (v.log(x)) • Then EXP. is called for 
computation of expfu.Iogfxl) . The result is calculated according 
to the formula and Is returned to the calling program. 


3!, ERROR ANALYSIS. 

40^1,001* pairs (x,zl of random numbers were generated with 
distribution the product of a right half of a Cauchy distribution, 
and the product of two Cauchy distributions. x**z was computed for 
each of these oairs, first using the routine, and then using double- 
precision operations. The maximum absolute value of the relative 
error in the routine was found to be 7.196 * 10**(-10» for these 
pairs. 


4. EFFECT OF ARGUMENT ERROR. 

If a small error e(xl occurs in the base x, and a small error efzl < 
= “'fxlt i.e*(y)) occurs in the exponent z, the error in the result 
w is given approximately by 

w * (log X * efzl ♦ z * efxl/xl . 
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1, ROUTINE'S FUNCTION 


1.1» Type. a FORTRAN exponentiating routine. It accepts an 
argument set consisting of a complex base and a fixed-point 
exponent* and returns a complex result. 

i.?. Purpose. To accept calls by name for ZTOI* generated by 
FORTRAN programs raising complex quantities to fixed-point 
exponents. 


?, method. 


See the description of ZTOI. for the algorithm. The argument set is 
checked upon entry. It is invalid if either argument is infinite or 
indefinite* or if the base is zero and the exponent is not greater 
than zero. In these cases* POS. INOEF. is returned and a 
diagnostic message is issued. Otherwise the result of the 
computation is returned to the calling program. 


T. ERROR ANALYSIS. 


Not applicable* since the only errors are round-off errors. 


4. effect of argument error. 


If a small error e* occurs in the base b, the error in the result is 
given approximately by n * b ♦♦ fn-1) * e** where n is the exponent. 
The absolute value of this expression is approximately the absolute 
error. If the error e* is significant* the error in the result 
should be found by substitution of the possible argument values in 
the expression b ** n. 
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ROUTINE t ZTOI. 

1. ROUTINE-S FUNCTION. 


1.1. Type. A FORTRAN exponentiating routine. It accepts an 
argument set comprising a complex and a fixed-point argument, 
and returns a complex result. 

1.?. Purpose. To accept calls by value for ZTOI. . generated by 
FORTRAN programs which raise complex quantities to fixed-point 
exponents. 


?, MFTHOO. 


Let b be the base and pC>OI the exponent. If p has a binary 
representation 000...01(n)ifn-l)...if1)i(0) where each i(J)C0<J<nl 
is 0 or 1, then 

P = i(0) . 2® * i(l». 21 *...*i(n) . 2 ** n 

and n = tlog(2Jp3 = greatest integer not exceeding log(2}p. Then 
b ** p = Prod Cb *» 2 ** I l 0<l<r 8 , i (J ) = 11 , 

The numbers b = b 2 *»o . b® ♦ b* «••.« b**2**n are generated by 
successive squarings, and the coefficients ifOI...., ifn) are 

obtained as sign bits of successive circular right shifts of p 
within the computer. A running product is formed during the 

computation, so that smaller powers of b may be discarded. Thus, 
the computation becomes an iteration of the algorithm 
b ♦* p = 1 if paO 

b ♦♦ p = (b*) ♦♦ p/2 if p>0 and p is even 
b p = b. (b*) ** (p-11/2 if piO and p is odd 

Upon entry, if the exponent p is negative, p is replaced by -p and a 

sign flag is set, b»*p is computed according to this algorithm, and 
if the sign flag was set, the result is reciprocated, before being 
returned to the calling program. 

The input range is the collection of pairs of bases b and exponents 
p such that b is non-zero if p is negative, both arguments are 
definite and in-range, and the result is in-range. 


3. ERROR ANALYSTS - not applicable. 


4. EFFECT OF ARGUHENT ERROR. 


If a smalI error e* occurs in the complex base b, the error in the 
result is given approximately by p ® b**<p-ll * e*• If e* is 
significant, the absolute value of the error in the result is less 
than or equal to 

Ipl » (Ibl ♦ lb + e*ll**fp-ll * le'l . 
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The Mathematical Library routines are classified according to the 
following criteria. 

f?outinei the name of the SCGPE deck concerned. 

Fntryl the name of the routine's entry ooint. 

Cal 1st by name or by value. 

Checkingt of arguments by the routine. 

Argumentt the type of arguments to the routine. 

FL = floating-point 
FI = fixed-point 
0 = double-precision 

C = complex 
A = any 


Pesulti 

the 

tyoe of 

the result 

(or results) 

• 


Funct ion I 

whether an 

external (Ext) or intrinsic (Int) 

function. 

Routine. 






EUDSliflDA 

A LOG 

ALOG 

Name 

Yes 

FL 

FL 

Ext 


ALOGIO 

Name 

Yes 

FL 

FL 

Ext 

ATAN 

ATAN 

Name 

Yes 

FL 

FL 

Ext 

ATAN? 

ATAN? 

Name 

Yes 

(FL.FL) 

FL 

Ext 


ATAN?. 

Val ue 

Yes 

(FL.FL) 

FL 

Ext 

ATANH 

ATANH 

Name 

Yes 

FL 

FL 

Ext 

CCOS 

CCOS 

Name 

Yes 

C 

C 

Ext 

CEXP 

CEXP 

Name 

Yes 

C 

C 

Ext 

CLOG 

CLOG 

Name 

Yes 

C 

C 

Ext 

COG 

COS 

Name 

Yes 

FL 

FL 

Ext 

COSO 

COSO 

Name 

Yes 

FL 

FL 

Ext 

COSH 

COSH 

Name 

Yes 

FL 

FL 

Ext 

CSIN 

CSIN 

Name 

Yes 

C 

C 

Ext 

CSOPT 

CSQRT 

Name 

Yes 

C 

C 

Ext 

OACOS 

OACOS 

Name 

Yes 

0 

0 

Ext 

OASTN 

OASIN 

Name 

Yes 

0 

0 

Ext 

OATAN 

OATAN 

Name 

Yes 

0 

0 

Ext 

OATAN2 

OATAN? 

Name 

Yes 

(0,0) 

0 

Ext 

OCOSH 

OCOSH 

Name 

Yes 

0 

0 

Ext 

OEXP 

OEXP 

Name 

Yes 

0 

0 

Ext 

OLOG 

OLOG 

Name 

Yes 

0 

0 

Ext 

OLOGIO 

OLOGin 

Name 

Yes 

0 

D 

Ext 

OMOO 

OMOO 

Name 

Yes 

(0,0) 

0 

Ext 

OHOO. 

OMOO. 

Va 1 ue 

No 

(0,0) 

0 

Ext 

OCOS 

OCOS 

Name 

Yes 

0 

0 

Ext 
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Sasstsinaj. 

Argument. 

Result. 

EanaiiaiiA 

OSIN 

OSIN 

Name 

Yes 

0 

0 

Ext 

DSINH 

DSINH 

Name 

Yes 

D 

0 

Ext 

nSQRT 

OSQPT 

Name 

Yes 

0 

D 

Ext 

OTAN 

OTAN 

Name 

Yes 

0 

0 

Ext 

OTANH 

OTANH 

Name 

Yes 

0 

0 

Ext 

OTOn* 

OTOO$ 

Name 

Yes 

C0,0) 

0 

• 

OTOI» 

DTOIt 

Name 

Yes 

(0,FI» 

0 

- 

OTOX* 

OTOX$ 

Name 

Yes 

(0,FLI 

0 


OTnz» 

OTOZf 

Name 

Yes 

(0,0 

c 

- 

FRF 

ERF 

Name 

Yes 

FL 

FL 

Ext 

FRFC 

ERFC 

Name 

Yes 

FL 

FL 

Ext 

FXP 

EXP 

Name 

Yes 

FL 

FL 

Ext 

TTOO* 

ITOOt 

Name 

Yes 

(FI,0) 

0 


TTOJ* 

ITOJf 

Name 

Yes 

(FI,FI> 

FI 


TT0X» 

ITOX$ 

Name 

Yes 

(FI,FL» 

FL 

- 

ITQZ» 

ITOZ* 

Name 

Yes 

(FI,0 

C 

- 

SINCOS. 

SIN 

Name 

Yes 

FL 

FL 

Ext 


COS 

Name 

Yes 

FL 

FL 

Ext 


SIN. 

Value 

Yes 

FL 

FL 

Ext 


COS. 

Va 1 ue 

Yes 

FL 

FL 

Ext 

SINO 

SIND 

Name 

Yes 

FL 

FL 

Ext 

SINH 

SINH 

Name 

Yes 

FL 

FL 

Ext 

SORT 

SORT 

Name 

Yes 

FL 

FL 

Ext 


SORT. 

Vat ue 

Yes 

FL 

FL 

Ext 

TAN 

TAN 

Name 

Yes 

FL 

FL 

Ext 

TANO 

TANO 

Name 

Yes 

FL 

FL 

Ext 

tank 

TANH 

Name 

Yes 

FL 

FL 

Ext 

XTOO* 

XTOO$ 

Name 

Yes 

(FL,0 

0 

- 

XTOI» 

XTOTt 

Name 

Yes 

(FL,FI) 

FL 

- 

XTOY* 

XTOV* 

Name 

Yes 

{FL,FL) 

FL 

- 

XTOZ* 

XTOZ$ 

Name 

Yes 

(FL,C» 

C 


ZTni» 

ZTOI$ 

Name 

Yes 

(C,FT» 

C 

- 

RAMF 

RANF 

Name 

No 

A 

FL 

Ext 


RANGET 


No 

FL 

- 

Subroutine 

ranset 

RANSET 

Name 

No 

FL 

mm 

Subroutine 

A NO 

AND 

Name 

No 

(A,A, ...1 

A 

Int 

COM PL 

COMPL 

Name 

No 

A 

A 

Int 

LOCF 

LOCF 

Name 

No 

A 

FI 

Int 

MASK 

MASK 

Name 

Yes 

FI 

A 

Int 

OR 

OR 

Name 

No 

(A,A,...) 

A 

Int 

SHIFT 

SHIFT 

Name 

No 

(A,FI» 

A 

Int 

XOR 

XOR 

Name 

No 

(A,A, ...I 

A 

Ext 

COUNT 

COUNT 

Name 

No 

A 

FI 

Int 

AOS 

ABS 

Name 

No 

FL 

FL 

Int 


TABS 



FI 

FT 
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Ron tine. 

£Qi£:L>. 


£!3££lsJUldA 

Argument. 

Result. 

FunctIon. 

AIMAG 

AIMAG 

Name 

No 

C 

FL 

Int 

flINT 

AINT 

Name 

No 

FL 

FL 

Int 

AMAXO 

AMAXO 

Name 

Nt> 

CFI,FI,...) 

FL 

Int 

AHAXl 

A MAXI 

Name 

No 

<FL,FL,...J 

FL 

Int 

AMINO 

AMINO 

Name 

No 

IFI.Fl,...1 

FL 

Int 

AMINt 

AMINl 

Name 

No 

CFL,FL,...) 

FL 

Int 

AMOO 

AHOO 

Name 

No 

<FL,FL) 

FL 

Int 

CMPLX 

CMPLX 

Name 

No 

«FL,FL» 

C 

Int 

CONJG 

CONJG 

Name 

No 

C 

C 

Int 

DABS 

DABS 

Name 

No 

0 

0 

Int 

OBLE 

09LE 

Name 

No 

FL 

0 

Int 

DIM 

DIM 

Name 

No 

CFL,FL) 

FL 

Int 

OMAXl 

OMAXl 

Name 

No 

(OfO« ...1 

0 

Int 

DMINl 

DMINl 

Name 

No 

fO. Of . . .) 

0 

Int 

OSIGN 

OSIGN 

Name 

No 

(0,0) 

0 

Int 

FLOAT 

FLOAT 

Name 

No 

FI 

FL 

Int 

lOIM 

lOIM 

Name 

No 

(FI,FI) 

FI 

Int 

I NT 

I NT 

IFIX 

IDINT 

Name 

No 

FL 

FI 

Int 

ISIGN 

ISIGN 

Name 

No 

(FI,FT) 

FI 

Int 


SIGN 



(Fl,FL) 

FL 


MAXO 

MAXO 

Name 

No 

(FI,FI, ... ) 

FI 

Int 

MAXI 

MAXI 

Name 

No 

(FL,FL,...) 

FI 

Int 

MINO 

MINO 

Name 

No 

(FI,FI,...) 

FI 

Int 

MINI 

MINI 

Name 

No 

(FL,FL,...) 

FI 

Int 

MOB 

MOD 

Name 

No 

(FI,FI) 

FI 

Int 

REAL 

REAL 

Name 

No 

C 

FL 

Int 


SNGL 



0 

FL 


A COSIN. 

ACOS 

Name 

Yes 

FL 

FL 

Ext 


ASIN 

Name 

Yes 

FL 

FL 

Ext 


ACOS. 

Value 

Yes 

FL 

FL 

Ext 


ASIN. 

Val ue 

Yes 

FL 

FL 

Ext 

ALOG. 

ALOG. 

Value 

Yes 

FL 

FL 

Ext 


ALOGID. 

Va 1 ue 

Yes 

FL 

FL 

Ext 

ATAN. 

ATAN. 

Value 

Yes 

FL 

FL 

Ext 

ATANH. 

ATANH. 

Va 1 ue 

Yes 

FL 

FL 

Ext 

CABS. 

CABS 

Name 

Yes 

C 

FL 

Ext 


CABS. 

Va 1 ue 

Yes 

c 

FL 

Ext 

ccos. 

CCOS. 

Value 

Yes 

C 

C 

Ext 

CEXP. 

CEXP. 

Value 

Yes 

C 

c 

Ext 

CLOG. 

CLOG. 

Value 

No 

c 

C 

Ext 

COS=STN 

COS.SIN 

Value 

No 

FL 

(FL,FL) 

Helper 

CSIN. 

CSIN. 

Value 

Yes 

C 

C 

Ext 

CSORT. 

CSQRT, 

Value 

No 

c 

C 

Ext 


£2yliaa.t 

£Dica.t 


£l3gCbiD3a 


Essyilx 

Function. 

OASNCS. 

DAGOS. 

Value 

Yes 

0 

0 

Ext 


DASIN. 

Value 

Yes 

0 

0 

Ext 

OATAN. 

DATAN. 

Value 

Yes 

0 

0 

Ext 

OATANZ. 

OATAN2. 

Value 

Yes 

f0,0) 

0 

Ext 

OEXP. 

DEXP. 

Value 

Yes 

0 

0 

Ext 

OHYP. 

DCOSH. 

Value 

Yes 

0 

0 

Ext 


OSINH. 

Value 

Yes 

0 

0 

Ext 

OLNLOG. 

OLOG. 

Value 

No 

0 

0 

Ext 


OLOGin. 

Value 

No 

0 

0 

Ext 

OSNCOS. 

DSTN. 

Value 

No 

0 

0 

Ext 


OCOS. 

Va lue 

No 

0 

0 

Ext 

OSaRT. 

OSQRT. 

Value 

Yes 

0 

0 

Ext 

OTAN. 

OTAN. 

Value 

Yes 

0 

0 

Ext 

OTANH. 

OTANH, 

Value 

Yes 

0 

0 

Ext 

nino. 

OTOO. 

Value 

No 

<0,01 

0 

- 

OTOI. 

OTOI. 

Value 

No 

(0,FI) 

0 

- 

OTOX. 

OTOX. 

Value 

No 

C0,FLI 

0 

- 

HTOZ. 

OTOZ. 

Value 

No 

10,0 

c 


FRF, 

EPF 

Value 

Yes 

FL 

FL 

Ext 


ERFC. 

Val ue 

Yes 

FL 

FL 

Ext 

EXP. 

EXP. 

Value 

Yes 

FL 

FL 

Ext 

HYP, 

COSH. 

Value 

Yes 

FL 

FL 

Ext 


SINH. 

Value 

Yes 

FL 

FL 

Ext 

HYPERS. 

HYPERS, 

Value 

No 

FL 

(FL,FL) 

Helper 

ITOO. 

TTOO. 

Value 

No 

CFT,0> 

0 

• 

ITOJ. 

ITOJ. 

Value 

Yes 

fFI,FI) 

FI 

- 

ITOX. 

ITOX. 

Value 

No 

(FI,FLI 

FL 

- 

TTOZ. 

ITOZ, 

Value 

No 

tFI,Cl 

C 

- 

FINCSO. 

COSO. 

Value 

Yes 

FL 

FL 

Ext 


SINO, 

Value 

Yes 

FL 

FL 

Ext 

TAN, 

TAN, 

Va lue 

Yes 

FL 

FL 

Ext 

TANS. 

TANO. 

Value 

Yes 

FL 

FL 

Ext 

TANH. 

TANH. 

Value 

Yes 

FL 

FL 

Ext 

XT09. 

XTOO, 

Value 

No 

CFL,D> 

0 


XTOI. 

XTOI. 

Value 

Yes 

CFL.FTI 

FL 

- 

XTOY. 

XTOY, 

Value 

No 

CFL,FL) 

FL 

- 

XT07. 

XTOZ. 

Value 

No 

(FL,CI 

C 


ZTOI. 

ZTOI. 

Value 

No 

fC,FI> 

C 
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A1f routines in the FORTRAN coamon. library checking arguments and 
issuing error messages allow for standard and non-standard error 
recovery* as described in the FORTRAN Extended Version 4 Reference 
Manual. Routine* the name of the loader deck concerned. The 
structure of these routines satisflest 

Word II VFD 42/*<routine*s na«e>, 18/< relative position of entry 
point> 

When executing under traceback mode* register AO holds the field 
length when in the main program* and the first word address of the 
parameter list in the previous call, otherwise. In normal execution 
each routine must save the contents of AO before using this register, 
and before calling any other routine. AO’s contents must be restored 
upon return to the calling routine. 

The symbols SYSARG. and SYSERR. are two entry points in the FORTRAN 
common library utility package FORSYS. . A call at SYSARG. with a 
••bad* argument li.e., negative* zero* infinite or indefinite! in XI 
will return with X2 holding the address of the text of an appropriate 
error message. A call to SYSERR. with an error number in XI and the 
address of a diagnostic message in X2 will result in the printing of 
the diagnostic message and a traceback listing, provided that the 
first two words of each routine are as above* the return Jump to 
SYSERR. is in the upper half of a word, and the lower 18 bits contains 
a pointer from word 1 to the return Jump. 

The seauence of events on executing math library routines which 
issue diagnostic messages isi 
(a.I Enter routine. 

<b.) Check arguments. If valid, compute result and return 

through entry point. ISome routines also check the result 
before return.I If invalid, go to Cc.l. 

Cc.l Enter contents of register AO in TEMPAO. and enter the FWA 
of the parameter list Cnow in At) into AO . 

Id.) Call SYSARG* to obtain the address of an error message in 
X? , if the argument is infinite or indefinite (or zero or 
negative)’ In this case, go to If.). 

(e.) Otherwise, enter the address of an appropriate error 
message directly into register X2 • 

l f. ) Enter the error number into XI . (See the FORTRAN Extended 

Reference Manual.) (Step (f.) may precede step Id.).) 

l g. ) Return Jump to SYSERR. to initiate error actions. (Lower 

part of RJ word = trace pointer.) If non-standard error 
recovery is specified through a previous call to SYSTEMC , 
transfer will return to the supplied recovery routine. If 
standard error was inhibited* the Job aborts. Otherwise, 
control will return to the calling routine* at step (h.). 
(h.l The appropriate indefinite or infinite quantity is entered 
into X6 , and the contents of AO are restored from 
TEMPAO. . 

(1.) Return through the entry point. 
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A list of error numbers and diagnostic messages is given in the 
FORTRAN Extended Reference Manual* 

As the first routines to be rewritten in a project to implement 
full error checking in all routines* some routines (listed in Appendix 
A) now detect errors and issue messages for all bad arguments passed 
to them. These routines call new routines SYS=AIO or SYS=1ST (at 
entry points SYSAIO. or SYSIST* * respectively} for error processing. 
The sequence of events on executing these routines is* 

(a.I Enter routine. 

(b.l Check arguments. If valid* compute result and return 
through entry point. (Some routines also check the result 
before return.) If invalid* go to (c.). 

(c.) Set 82 with pointers indicating error number* partial 
message, and register residence of bad argument. The 

format is given in the method description of routine 
SYS=1ST . The partial message will be ignored if the 
argument is infinite or indefinite. 

(d.) Set up the arguments in registers XI , )(2 , X3 and X4 (or 
Just XI * X2 if one argument) according to the rules in 
section III of the Introduction. 

(e.l Return Jump to SYSIST, or SYSAIO. to initiate error 

processing, SYSAIO must be chosen if there is more than 
one argument. The return Jump must be in the upper 31? bits 
of a word. The next 12 bits are zero* and the next 18 bits 
must include a pointer to a trace word, as described above, 
(f.) Testing commences. A parameter list is built up from 

values in XI * X2 * X3 , X4 to allow non-standard error 
recovery. If the routine calling the routine calling 

SYS=AIO made this call in the format 
+ RJ =X<routine> 

- VFD 30/1 

go to step g below. Otherwise* set AO to point to the 
reconstructed parameter list* set XI to the error number* 
set X2 to the first word address of the constructed 

message, then execute the communication cell SYSAIO. * 
after traceback linkage information has been inserted in 
its lower 18 bits. 

(g.) Return +IN0. in registers X6 and X7 * and restore registers 
AO , XI , X2 (and X3 and X4 , if entry was to SYS=AID). 
s = .746926109335419 * 10»*-3 
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The times listed below were determined empiricalty* and arguments to 
routines were chosen as many as practicable to cover all the possibilities 
for times to each routine* CYBER 76 times were obtained through the machine 
instruction 016|0 which accesses a hardware clock* while CYBER 72* 73 and 74 
times were obtained by observing variations in speed of two equivalent loops 
in central memory* one of which catted the routine being timed* These 
variations in speed were obtained through use of a system-maintained real¬ 
time clock which is synchronized with a hardware clock on one of the data 
channels. These times do not include time for setting up arguments and 
parameter lists* but measure from the time a return-lump to the routine Is 
issued* to the time that the next instruction in sequence is issued* All 
times given are in minor cycles (or clock-periodsl* On CYBER 72* 73 and 74* 
1 minor cycle = 100 nanoseconds* while on CYBER 76* 1 clock-period = 27*5 
nanoseconds* On CYBER 171* 172 and 173* 1 minor cycle = 50 nanoseconds* 

Certain facts should be noted* On CYBER 76* a return lump may be delayed 
in execution if the instruction stack control has requested one or more 
instruction words that have not arrived at the instruction stack* Thus* 
CYBER 76 routine times depend on how the routine is called* On CYBER 72 and 
73* a floating instruction executes at least 4^ minor cycles faster if either 
of the operands is zero* infinite or indefinite* If in the course of 
evaluating an algorithm for computation of a function* a routine happens to 
produce an intermediate zero result* it will execute faster by at least 48 
minor cycles if this intermediate result is combined arithmetically with 
anything else* The number of possibilities for this case is too large for 
enumeration in this appendix* 

Some routines will naturally call others* but the time listed under each 
routine only the time spent in that routine * and does not include time spent 
in return Jumps to and execution of other routines* To find total execution 
time in a routine* one must add times for execution at entry points with 
arouments I isted after an “&*** 

Timings are supplied here for valid argument sets only* Take the time for 
the first alternative listed which covers the argument concerned* 


I 
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Routine 

Entry Points 

Arguments 

RTimes at Entry Points (argument) Times for CYBER 



173 

72 

73 

74 

76 

A8S 






A3S (Any valid) 


100 

79 

58 

66 

ACOSIN. 






ACOS (X) 






X va1id 



56 

35 

47 

% ACOS. (x) 






ASIN (x) 






X valid 



59 

35 

41 

t ASIN. (xl 






ACOS. (xl 






X valid and! 

X s 0. 

812 


741 

159 

119 

X = 1. 

306 


234 

127 

85 

X * -1 • 

307 


237 

127 

87 

X in (-.5*.5) 

05 0 


897 

159 

116 

X not in (-.5..5)* time 
= a*b*n where n is 
the loop count, as defined 
in the ACOSIN. description. 






a = 



1138 

207 

147 

b = 



114 

18 

12 

X in (“I..,*.?), 
add to |xj timet 



10 

5 

4 

ASIN. (X) 






X valid and! 

X = 0 • 

823 


763 

153 

123 

X = 1. 

292 


220 

120 

87 

X = -1 • 

293 


219 

120 

90 

X in (-.5,.5) 

|x| in (.S,l.), time 
= a*-b*n, where n is 
defined in the ACOSIN* 
description. 

958 


904 

152 

126 

a = 



1170 

226 

168 

b = 



115 

15 

12 

AIMAG 






AIMAG (Any valid) 


101 

81 

54 

62 

AINT 






AINT (Any va! id) 


121 

98 

66 

6 0 
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Routine 

Entry Points 

Arguments 

ITiwes at Entry Points largunentl Times for CYBER 

173 72 73 74 76 


ALOG 


ALOGIB lx» 



67 

42 

46 

1 ALOGIO. lx) 






ALOG lx) 



67 

40 

49 

t ALOG. lx) 






ALOGIO. lx) 






X Infinite or indefinite 



253 

192 

110 

1 SYSAIO. CApoend. 

8) 





8. 



266 

199 

120 

% SYSAIO. CAppend. 

B) 





X vafid«x<0. 



285 

213 

129 

1 SYSAIO. CAppend. 

B) 





X valid. X a y*2**n* 






n integral. l<v<2. and 






l<y<1.1072 


860 

892 

179 

129 

1.1072<y<1.3572 


860 

892 

176 

129 

1.3572<y<1.6072 


861 

891 

177 

128 

1.6l!72<y< 1.8572 


860 

892 

179 

129 

1.8572<v<2 


990 

1012 

212 

141 

ALOG. (x) 






X infinite or indefinite 



298 

176 

97 

t SYSAIO. CAppend. 

B) 





0. 



311 

183 

112 

1 SYSAIO. CAppend. 

B) 





X valid.x>0 



330 

192 

117 

i SYSAIO. CAppend. 

B) 





X valid.xa y*2**n. n integral 





t.<y<1.8572 


941 

814 

197 

119 

1.857?<2 


1072 

933 

218 

143 
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Routine 

Entry Points 

Arguments 

&Times at Entry Points (argument! Times for CYBEP 

173 7? 73 74 76 


AMAXQ 


AMAXQ (x(l)»..., x(n!l 
n-2 

240 

178 

121 

112 

n=3 

338 

250 

159 

140 

each add* 

99 

73 

42 

32 

AMAXl 

AMAXl (x(ll*..*, x(n!) 

n=2 

232 

178 

104 

106 

each add* 

110 

83 

45 

14 

AMINfl 

AMINO ((x(l)! 

***« x(n) 
n=2 

237 

179 

112 

105 

n=3 

338 

252 

148 

113 

each add* 

100 

72 

43 

32 

AMINl 

AMINl (x(ll«**.* x(nl) 

227 

179 

108 

105 

n = 3 

333 

252 

163 

142 

n=4 

436 

328 

189 

172 

each add* 

105 

78 

44 

34 

AMOn 

AMOO (x,v! 
y#0 

248 

207 

111 

133 

ANO 

AND (x(l),..** x(n)} 

n=2 

217 

163 

103 

103 

n=3 

282 

212 

112 

118 

n=4 

347 

262 

133 

141 

each add. 

65 

49 

22 

19 
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Routine 

Entry Points 

ftrguments 

S,Tlmes at Entry Points 


flTAN 

ATAN fxl 

^ ATAN. (xl 

ATAN. (xJ 

X valId Ix|<l. 

X valid |x|>l. 

ATAN? 

ATAN2 <y,xl 

4 ATAN2. (xl 

ATAN2. (y,xJ 

(y«xl valid and ... 

x=0»v#0. 

x#O,V=0 

lx|>|vl>0 

ly!>|x!>0 

ATANH. 

ATANH.fxl 

X valid andt 
x»0 

.75<x<1.5 

x>1.5 

CABS. 

CABS (zl 

7 valid 

I CA3S. (zJ 

CABS. (x*-i»yl 

x+i*y valid 
and x=y=0. 
x#0. or v#0. « 
special case. (See 
routine’s descriotion) 
and otherwise valid 

CCOS 

CCOS (zl 

z valid 

4 HYPERB. (Ifi(z) 

4 COS.SIN (re(zll 
|i«(zl|> 741.67 

4 SYSERR. (Append. Bl 


(argument 1 

173 

Times 

72 

for CYBER 

73 74 

76 



66 

32 

53 

1059 


756 

187 

141 

1092 


784 

203 

201 



78 

53 

78 

898 


850 

246 

190 

981 


835 

2 76 

187 

1167 


1085 

249 

161 

1185 


1077 

241 

172 


682 

908 



203 

202 




105 

38 

43 

276 


225 

138 

85 

715 

715 


786 

664 

283 

283 

197 

181 


546 

436 

180 

119 


468 

348 

363 

131 
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Routine 

Entry Points 

Arguments 


SiTimes at Entry Poirtts largunent) 

Times 
173 72 

for 1 
73 

CYBER 

74 

76 

cnos. 

CCOS. fzl 

z valid 

327 

279 

78 

71 

t HYPERB. (infzM 

8. COS.SIN freCzl ) 

CEXP 

CEXP (zl 

1refz»l> 741.67 

356 

273 

158 

120 

1 SYSERP. fAppend. B) 
z valid 

487 

403 

155 

115 

1 EXP, (re(z)l & COS.SIN (ImlzM 

CEXP. 

CEXP. (zl 

z valid 

262 

225 

74 

60 

1 EXP. (relzl) S. COS.SIN (imlzll 

CLOG 

CLOG (zl 

Z=0. 

291 

213 

IIS 

76 

L SYSARG= SYSERR. (Append. B) 
z valid 

163 

131 

102 

69 

& CLOG, (zl 

CLOG. 

CLOG. (zl 

z valid 

253 

199 

95 

50 

1 ATAN2, ((ifR(z), re(zll} 
i CABS, (zl t ALOG. (Izll 

CMPLX 

CMPLX (x,yl 

x.y valid 

126 

103 

64 

84 

COM PL 

COMPL (xl 

03 

69 

55 

54 
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Routine 

Fntry Points 

arguments 

STimes at Entry Points Cargumentl Times for CYBER 

17T 72 73 74 76 


CONJG 

CONJG (z) 

z vaI i d 


COS. — see SINCOS. 


128 101 S8 68 


COSH. -- see HYP. 


C0S = SIN 

COS.SIN fx) 


1x1 > Di*2**® 


307 

244 

108 

90 

ix|=v(fflOd2i3i) , 

0<y<2Di 0<v<pi^4 

1463 

1561 

1380 

242 

215 

Di/4<y<pi/2 

1715 

1880 

1649 

269 

234 

Di /2<v53pi/4 

1716 

1879 

1649 

269 

2 34 

3p i /4< v<p i 

1734 

18 85 

1655 

282 

245 

oi<y?5pi/4 


1884 

1657 

323 

245 

SDi/4<v<3pi/? 


1886 

1659 

319 

244 

3pU2<y<7pi74 


1887 

1658 

319 

244 

7p i74<v<2Pi 

1693 

1885 

1635 

267 

232 
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Routine 

Entry Points 

Arguments 

&Times at Entry Points (argument) Times for CYBER 

173 72 73 74 76 


COUNT 





COUNT (x) 

148 

133 

49 

62 

CEIN 

CSIN (zl 

|re(z) 1 > Di*2*»« 

366 

295 

162 

121 

1 COS.SIN (re(z)) 

1 SYSERR, (Apoend, B) 
lim(z)1 > 741.67 

470 

362 

221 

136 

1 COS.SIN (re(z)) 

1 SYSERR. (Aooend. B) 
z valid 

551 

436 

181 

123 

8, COS.SIN (re(z)) 

1 HYPERB. Cimlzl) 

CSIN, 

CSIN. (z) 

z valid 

315 

248 

77 

79 

8 COS.SIN (relz)) 

8 HYPERB. (im<z)l 

CSORT 

CSQRT (z) 

z valid 

153 

115 

93 

67 

8 CSQRT. (z) 

CSORT, 

CSQRT, (z) 
z=0. 

287 

219 

103 

58 

8 CABS. (0.) 

8 SORT, 10.) 
z valid* zi^O 

477 

376 

265 

90 

8 CABS, (z) 

8 SQRT. (l/2(lzl ♦ Ire(z)t)) 

DABS 

0A8S (X) 

X valid 

144 

111 

70 

72 
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Routine 


Entry Points 

Arguments 

^Times at Entry Points largumentl Times for CYBER 

173 72 73 74 


OASNCS. 

HACOS. 

0<x<.09375 
.09375<x<.7071 
.701<x<.9956 
.9956<x<l 

DASTN. 

0<x<.09375 
.09375<x<.7Q71 
.701<x<.9956 
.9956<x<l 


3344 

529 

4844 

853 

4823 

841 

4228 

756 


3260 

492 

4756 

814 

4779 

820 

4197 

7 36 


76 
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Routine 

Entry Points 

arguments 

STimes at Entry Points (argument} Times 

173 12 


for CYBER 
73 74 



OAT AN 

DATAN (x> 

X valid 130 42 143 

1 OATAN. (X) 

DATAN. 

OATAN. 

X valid* and: 

|xl < 1. 144 74 40 

1 OTN. (see routine* description) 

|x| > 1. 320 134 73 

1 DATCOM. (see routine* description) 

0ATAN2 

OATAN? fy.x) 

y.x valid, and (y,x)E(Q,0) 124 46 66 

1 OATAN2. ((y.x)) 

OATAN2. 

OATAN? (y.x) 

where both are valid, and 
(y.x)^(O.O). and: 

|y| < |x| 276 144 65 

1 OATCOM. (see routine* description) 

|y| > lx I 283 175 71 

& OATCOM. (see routine* description) 

OATCOM. 

OATCOM. (y.x) (from 0ATAN2. ) 

argument set validated. If n 
is nearest integer to 
8*min(|x|.1y|)/max(|x| .lyl), 
then: 


n = 0 

315 0 

5 21 

337 

nEO and min( I x | » ! y |)-n/8»max (| x | , | y I) EO 

3735 

664 

417 

otherwise 

3725 

663 

417 


OTN. 

y (from OATAN. ). valid. 

Tf n is nearest integer to 
8»y, then: 


r>=Q 

2736 

451 

287 

nEO and (y » n/8)E?) 

3356 

5 87 

367 

otherwise 

1212 

307 

200 
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Rou tine 

Entry Points 

Arguments 

ITimes at Entry Points largumentl Times 

173 72 


for CYBER 
73 74 


76 


08LE 

OBLE <xl 


X valid 

□COS 

OCOS (xJ 

98 

78 

52 

54 

X valid 

S DCOS. Cx) 

OCOSM 

□COSH (X» 

144 

121 

71 

67 

X valid 

i OCOSH. Ix» 

OEULER. 

□EULER, 


130 

52 

45 

(See description of routine OEULER* > 

DEXP 

□EXP (x) 


3719 

623 

361 

X vaiid 

1 DEXP.(x) 

□EXP. 

DEXP. (x» 

X valid andt 


117 

45 

49 

x<-643.2405fi355962g247139191409t 

K OEULER. 


515 

163 

107 

X otherwise* 


378 

147 

100 


S. OEULER. 
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Routine 

Entry Points 

Arguments 

STimes at Entry Points CargumentJ Times for CYBER 

173 72 73 74 76 


OHYP 

OCOSH. fx) 

X valid andi 


absCxI>42,3606303797914263858556020791 

560 

215 

127 

1 DEULER. 
absCxI/log 2 > 

48* 

546 

182 

101 

1 DEULER. 
absCxI/log 2 > 

24* 

658 

203 

125 

& DEULER. 
x in [-1/2 log 

2*1/2 log 2]* 

233 

136 

86 

1 DEULER. 

X otherwise 


719 

229 

132 


& DEULER. 


OSTNH. Cxi 

X vatid and* 


absCx)>42,360630379701426385855602079* 

STS 

202 

119 

1 DEULER. 
absCxI /log 2 * 

48* 

SIS 

160 

93 

t DEULER. 
absCxI/log 2 > 

24* 

625 

206 

136 

1 DEULER. 
X in [-1/2 log 

2,1/2 log 21* 

155 

94 

64 

1 DEULER. 

X otherwise 


720 

226 

134 


1 DEULER. 


DIM 

DIM fx,yl 

xty valid 191 150 84 96 
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Routine 

Entry Points 

Arguments 

STimes at Entry Points (argument) Times for CYBER 

173 72 73 74 76 


DNLOG. 

OLOGHl. (x) 

x=<2**n)*y 



l/2<y< 


7104 

7931 

6946 

1220 

7 61 


l/2»*.5<y<l 


6962 

7799 

6802 

1221 

762 

OLOG, 

(X) 

x=(2**n)*y 
l/2<y< 1/2**.5 


6797 

7576 

6631 

1158 

731 


l/2**.5<v<l 


6636 

7444 

6487 

1144 

731 

DLOG 








□ LOG 

(X) 

x=13. 



284 

215 

136 

8 3 


1 SYSARG. SYSERR. (Append. 

B) 







x< 0 



332 

251 

142 

105 


& SYSARG. SYSERR, (Append. 

B) 







X valid 



150 

96 

101 

68 


t OLOG. (x) 







OLOGIO 








OLOGlfl (x) 








x-0. 



284 

216 

130 

89 


t SYSARG. SYSERR. (Append. 

B) 







X* x<0 



333 

255 

216 

105 


1 SYSARG. SYSERR. (Append. 

B) 







X valid 



177 

144 

97 

68 


t OLOGIO. (x) 







OMAXl 








□MAXI 

(x(l),x(2)) 



909 

675 

320 

135 

OMINl 








□ MINI 

(x(i).x(2)) 



863 

644 

310 

134 

OMOO 








OMOO 

(x,yl 

X valid. v=0 



332 

243 

137 

77 


1 SYSARG. SYSERR. (Append, 

B) 







(x.y) valid 



266 

203 

34 

97 


I OMOO. (x,y) 

OMon. 

OMOO, <x,yl 


x.y valid.y#0 lx/vl>2«« 

2007 

582 

tx/yl??** 

1426 

431 

1 x/yl <?**« 

841 

281 


174 
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Routine 

Entry Points 

Arguments 

ITimes at Entry Points (argument) Times for CYBER 

173 72 73 74 76 


OSIGN 

OGIGN (x,y) 

X, V any 

OSIN 

DSTN (x» 

X valid 

S. OSIN* (x) 

OSINH 

OSINH Cx) 

X valid 

1 OSINH. (x» 


205 157 81 101 

162 102 83 76 

124 52 43 


OSNCOS. 

OCOS. (x) 


|x|>pl.2«* 


605 

501 

181 

129 

x=y(mod2pi)* 0<y<2ol 

0<y<pi/4* 

4671 

5129 

4475 

778 

516 

Pi/4<y<pi/2* 

5140 

5703 

4971 

8 44 

563 

oi/2<v<3pi74. 

5140 

5703 

4971 

846 

563 

3pi74<y<pi* 

5059 

5679 

4904 

851 

558 

pi<y<5pi74, 


5658 

4923 

920 

558 

5pi/4<y< 3pi72, 


5703 

4980 

908 

563 

3pi/2<y< 7pi/4 


5722 

4971 

909 

563 

7pi74<y<?pi 

5063 

5677 

4904 

850 

563 

OSIN. (X) 

|x|>Di .2»^ 


624 

511 

181 

137 

x=v(fflOd2pi). 0<y<2pi, 

0< y<pi/4, 

4750 

5093 

4446 

786 

520 

pi/4<v<oi7?, 

5078 

5695 

4933 

867 

566 

pi72<y<3pi74. 

5083 

5689 

4904 

864 

571 

3pi/4<y<p iy 

5139 

5715 

4971 

856 

575 

PiSy<5pi74. 


5715 

4980 

924 

571 

5pi74<y< 3pi72, 


5689 

4933 

934 

566 

3Di72<v< 7oi74 


5687 

4904 

935 

566 

7oi74<y<2Dl 

5141 

5718 

4980 

853 

571 
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Routine 

Entry Points 

Arguments 


ITimes at 

Entry Points (arquwent) 

Times 
173 72 

for 

73 

CYBER 

74 

76 

nSORT 

DSqPT (xl 
x<0. 


282 

234 

125 

85 

1 SYSARG. 
X valid 

SYSERR. (Append. 8) 

140 

107 

93 

60 


S OSQRT. <xl 

nSQRT, 

OSQRT. tx) 
x-Q* 

x=y*2**n 
n odd 
n even 
OTAN, 

HTAN. 

X valid andt 
x=0 

oi/4<x<oi/4 
Di /4<x<3pi/«» 

3pi/4<x<5pi/4» etc. 

5pi/*»<x<7oi/<»« etc. 

OTANH 

OTANH (x» 

X valid 12h 120 42 

t OTANH. (xl 


745 

746 


2371 

2247 

3663 

3474 

3666 


228 

231 


579 
579 
639 
6 33 
638 


HTA.NH. 

OTANH. ix> 


X valid andi 
lxl<l/8* 

i OEULER. (x) 

765 

217 

134 

!x|> 321 

214 

103 

62 

If X (or 2xJ=v*n*log(2)* n>47i 

1 OEULER. (2x1 

619 

163 

122 

otherwisei 

X OEULER. 12x1 

1055 

311 

171 


176 
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Routine 

Entry Points 

Arguments 



&Times at Entry Points (argument) 

Times 
173 72 

for 

73 

CYBER 

74 

76 

OTOO* 

DTOOS (x,yl 
(Q.tO. ) 

441 

341 

192 

158 


1. SYSERR. (Append. B» 

(0«y), to y>0 

3S2 

387 

153 

208 


(0»v)» to y<0 

439 

340 

195 

152 


8, SYSERR. (Append. B) 

x< 0 

410 

318 

169 

130 


8 SYSERR. (Append. B) 

(Xfy) valid 

863 

740 

236 

138 

OTon. 

OTOO. 

8 OLOG. (X) 

8 OEXP. (y*log xl 

(x,y) 

(0«y)« y>0 

114 

63 

66 

62 


x>0» x,y valid 

517 

466 

113 

79 

% OLOG. lx) 

& DEXP. Cy*log x) 

OTOI* 

0T0I$ (x,nl 
(a.fOl 

404 

312 

189 

136 


8 SYSERR. 

(?'.*nl.n<0 

418 

311 

195 

142 


8 SYSERR. (Append. B) 

(0.*n),n>0 

230 

188 

123 

192 


x> 0 

264 

231 

111) 

69 

oTor. 

OTOI. 

8 OTOI. (x,n) 

(x.n) 

if n<0« add* and replace 

467 

415 

114 

73 


n with -n 
(x,OI 

83 

65 

51 

51 


(x,l) 

364 

301 

126 

90 


(x*2» 

672 

575 

190 

128 


if n>2*timest* a(l)tb(ll< 
logl?)n<t<a(2» /t-bt?) log(2»n 
a(l)s 

380 

316 

227.9 

94.3 


a(2) = 

69.3 

14.5 

111.6 

105 


b(l) = 

292. 

257. 

38.8 

24.2 


b(2) = 

530 

489 

41.9 

33.6 
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Rou tine 

Entry Points 

Arguments 

S,Times at Entry Points largument) Times for CYBER 

173 72 73 74 


76 


OTOX* 

DTOX55 (x,y» 

( 0 *^ 0 ) 

1 SYSERR. CAppend, Bl 
(0«»yl♦y<0 

(O.iyJ tv<0 

1 SYSERR. 

x<0 

& SYSERR. (Append. B) 


(x.y) valid. 

x>0 

A OLOG. 

(X) 

A DEXP, 

(y*loq 2) 


DTOX. 

DTOX. (x,y) 
(O.ty) 

(x.y) valid 
I OLOG. 
i OEXP, 


(xl 

(y»log 21 


OTOZ* 


(x* 

z) 




(n.. 

0. 

♦i.O.I 




A 

SYSERR. 

(Append. 

81 

x<0 






A 

SYSERR, 

(Append. 

B) 

(0., 

zl 

, Pe(z)>0 




A 

SYSERR. 

(Append. 

B) 

(Q., 

zl 

, re(z)<0 

im(zl#0 



A 

SYSERR. 

(Append . 

B) 

(0., 

z) 

* re<z)<0 

im(zJsO 



A 

SYSERR. 

(Append. 

B> 

(x,z 

1 

valid 




1 ALOG. (X) 
t EXP. Ire(zl*l og xl 
% COS.SIN (im(z}«t 09 x1 


0T0 7. 

0T07. (x,Zl 
x»0. 

x«z valid* x^O 

A ALOG. (xl 
A EXP. (r«(zf»lo 9 x» 

A COS.SIN (im(z)«lo 9 xl 


426 

337 

199 

184 

2B9 

239 

245 

228 

4 26 

335 

107 

178 

383 

299 

176 

145 

708 

606 

236 

158 


95 

74 

59 

54 

460 

415 

85 

63 


403 

311 

189 

148 

342 

263 

168 

117 

277 

211 

102 

136 

432 

312 

221 

136 

432 

312 

223 

136 

762 

636 

231 

90 


103 

81 

63 

59 

480 

426 

149 

85 


178 
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EPF. 


EPF. (x) 


x<-E.F25 or -inf 

526 

189 

-F.625<x<-.477 

3094 

489 

-.477<x<0 

1172 

2 34 

x=0 

904 

230 

’1<x<.477 

1172 

235 

.477<x<5.625 

3090 

495 

x>5*635 or +101 

ERFH. {xJ 

527 

185 

x<-5.F?5 or -inf 

588 

213 

-5.62F<x<-.477 

3155 

518 

-.477<x<0 

1234 

255 

x=0 

965 

252 

?l<x<.477 

1234 

253 

.477<x<8 

3154 

513 


x> 8 

X infinite 
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Routine 

Entry Points 

Arguments 

&Times at Entry Points 


EXP 

EXP CxJ 

& EXP. Cx) 

EXP. (x> 

X infinite 

8, SYSAIO. CAppend. 8J 
X indefinite 

8 SYSAIO. (Append. B) 
X valid,x>741.67 

8 SYSAIO. (Append. Bl 
x valid x>512. 

X val id,x<-675.84 

8 SYSAIO. (Append. B) 
X valid,x<-51? 

X valid 

FLOAT 

FLOAT (x) 

X valid 

HYP, 

COSH, (xl 

X valid 
|x|<l/2 log 2 
X otherwise valid 

SINH. (x) 

X valid 
|x|<l/2 log 2 
X otherwise valid 

COSH (xl 

X valid 

STNH (xJ 

X valid 

HYPERB. 

HYPER8. (xl 

X valid, |x|<.22 
X valid, |xl>,22 
8 EXP. (xl 

TOIM 

lOIM (x,yl 

(x,yl valid 


(argument) 

173 

Times 

72 

for 

73 

CYBER 

74 

76 



34 

57 

38 



268 

140 

89 



201 

103 

58 



304 

155 

97 

932 


864 

184 

130 



298 

157 

119 

931 


865 

182 

140 

843 


804 

145 

112 


102 

82 

65 

56 


1296 

1385 


1313 

1426 

233 

233 

164 

167 

1325 

1457 


1351 

1498 

250 

257 

178 

177 



1495 

283 

200 



1559 

306 

214 

1649 

1772 

398 

1540 

311 

347 

136 

265 

95 


163 

127 

85 

103 
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Routine 

Entry Points 

Arguments 

KTimes at Entry Points (argument) Times for CYBER 

173 72 73 74 76 


I NT 

IFIX fx) 


X valid 

TNT 

lOINT 

ISIGN 


101 

81 

59 

56 

I^IGN (x,y» 

TTOO» 

ITODS (n,x» 


161 

125 

75 

96 

(fl • t 0 • ) 

1 SYSERR. (Append* 

8) 

365 

337 

164 

132 

(0*xl*x<0 

t SYSERR* (Append* 

Bl 

582 

483 

166 

123 

(0»x)*x>0 


496 

451 

109 

173 

n< 0 

S, SYSERR. (Append* 

6) 

322 

153 

128 

92 

n>0* x*log n overflows 

8 SYSERR* (Append* 

8 OLOG* (n» 

B» 

695 

584 

238 

914 

(n,xl validf n>0 

8 OLOG* (nl 

8 OEXP* Cx*log nj 

ITOO. 

ITOO* (n,x) 


598 

613 

264 

125 

(O.tx) 


146 

110 

78 

64 

(n*x) valid, n>0 

8 OLOG. (n) 

8 OEXP* (x»1og n) 


457 

397 

98 

80 
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Rou tine 

Entry Points 

Arguments 

STimes at Entry Points (argument) Times for CYBER 

17? 72 73 74 


ITOJ* 

ITDJS (»,n) 

1 ITOJ* (m^n) 


ITOJ. 

ITOJ* (m*n) 

m**n <2‘»« 


(mi 

*9} 

valid 

181 

95 

(mi 

► l)»m 

val id 

218 

131 

(«i 

1 ?) * m 

val i d 

283 

139 


if n>2,m>l* look at n in binary! 
for each 1 bit* add 
for each 0 bit, add 


76 


182 
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Pou tine 

Entry Points 

Arguments 

STimes at Entry Points (argumentl Times for CYBER 

173 72 73 74 76 


ITOX* 

ITOX$ (n,xl 


(OtO.I 


389 

267 

175 

158 

t SYSERR. (Append. 
(O.xl♦x>0 

B) 

352 

313 

114 

298 

(0.xl.x<0 


346 

268 

178 

149 

8, SYSERR. (Append. 

n< 0 

Bl 

289 

223 

135 

102 

8. SYSERR. (Append. 
n>0, 1x»lognI>741.67 

Bl 

459 

354 

246 

122 

8, ALOG. (nl 

S> SYSERR. (Append, 
(n.xl valid 

Bl 

315 

237 

245 

95 


«, ALOG. (n) 
t EXP. <x*log n) 

ITOX. 

ITOX. <n,xl 

(0,xJ 113 85 68 62 

(n.xl valid n>0 215 185 64 

I ALOG. (nl 
(n.zJ valid 

R ALOG. (n) 113 85 175 

& EXP. (x*1og n» 

ITOZ* 

TTOZJ (n,zl 


(0.0.4^i0. } 

1 SYSERR. (Append. 

B) 

376 

291 

165 

129 

(9»z)» re(z) <0tim(zl=(! 

& SYSERR. (Append. 

Bl 

395 

287 

199 

120 

(O.zl.re(z)>0 


238 

187 

210 

91 

im(zl=0 ( 0. zl . im (zl #0 . 

8. SYSERR. (Append. 

81 

376 

291 

165 

120 

re(zl<0 ln,zl.n<0 

& SYSERR. (Append. 

61 

316 

241 

144 

104 

(n.z) va!id 


632 

515 

211 

139 


\ ALOG. (nl 

& COS.SIN (im(z)*log nl 
& EXP. (re(z)*1og nl 

ITOZ. 

ITOX. (n,zl 84 42 24 

& XTOZ. (n,z) 
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Routine 

Entry Points 

Arguments 



ilimes at Entry Points (argument) 

Times 
1T3 72 

for f 
73 

CYBER 

74 

76 

LOCF 

LOCF 

(X) 

72 

60 

46 

49 

MAX!? 

?<AX0 

(xd)* ••• t x(n)) 

n=2 

222 

168 

105 

113 


n=3 

324 

240 

148 

134 


n=4 

422 

314 

191 

166 


each additional argument 

100 

73 

43 

31 

MAXI 

MAXI 

(xd) . ••• x(n) ) 

n=2 

249 

187 

111 

111 


r=3 

357 

270 

157 

141 


n=4 

467 

355 

202 

175 


each additional argument 

110 

83 

45 

34 

MASK 

MASK 

In) 

n>60 

263 

207 

111 

91 


8. SYSERR. CApoend. 8) 

n<0 

274 

210 

127 

83 


i SYSERR. (Append. 8) 
n valid 

181 

133 

103 

87 

MTNll 

MINO 

xtn)) 

n=2 

228 

169 

105 

102 


n=3 

328 

241 

148 

130 


ns:4 

429 

312 

191 

162 


each additional argument 

100 

72 

43 

28 

MINI 

MINI 

(x (1)f« x(n)) 

ns2 

242 

182 

110 

107 


n=3 

347 

259 

155 

137 


n=4 

454 

337 

199 

171 


each additional argument 

105 

77 

44 

38 

MOD 

MOO 

(x,y) 

(x.y) valid 

316 

268 

114 

133 
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PoutIne 

Entry Points 

Arquments 

iTimes at Entry Points 

(argument) 

173 

Times 

72 

for CYBER 

73 74 

?6 

OR 







OP xCnIl 







n=2 



209 

161 

103 

88 




274 

210 

124 

106 

n=4 



335 

258 

145 

130 

each additional argument 



63 

48 

21 

20 

PANE 







PANE {anything) 



189 

165 

63 

80 

RANGET (x) 



96 

79 

67 

66 

X will be modi fled 







PANSET 







RANSET (x) 



176 

134 

89 

82 

REAL 







REAL (u) 







SNGL (u) 



83 

69 

55 

54 

u valid 







SHIET 







SHIET (u,n) 







n valid 



128 

104 

60 

86 

SINCOS, 







SIN «x) 




64 

34 

42 

t SIN. (x) 







COS Cx) 




64 

34 

42 

1 COS. lx) 







SIN. <x) 







X Infinite or Indefinite 




169 

115 

75 

% SVSAID. (Append. 

B) 






x= 0. 


888 


821 

193 

141 

X validf |x|>pi*2**® 




166 

109 

79 

% SYSAIO. (Append. 

B) 






X valid. IxlEpi^?**® 


1283 


1256 

194 

141 

COS. (x) 







X infinite or indefinite 




169 

115 

75 

%. SYSAIO. (Append. 

8) 






0. 


831 


757 

188 

165 

X val id 1 X 1 Sp i »?‘»® 


1190 


1230 

188 

178 

X valid |x|>pi*2*** 





220 

165 

S, SYSAIO. (Append. 

B) 
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Routine 

Entry Points 

Arguments 

STimes at Entry Points (argument} Times for CYBER 

173 72 73 74 76 


SHRT 


SORT 

(X) 


78 

40 

37 


1 SORT, (xl 





SQRT 

. (x) 






X infinite, indefinite or negative 


222 

180 

279 


1 SYSAIO* (Apoend. B) 






X valid.x^O 

527 

523 

119 

101 


Q. 

244 

393 

196 

97 


SYEsAIO 
^ SYSAIO. 

(1 in lower half of RJ word! 359 133 

1 SYSERR. (Apoend* 3} 

(other than 1 in lower half of RJ word) 986 423 267 

K SYSERR. (Apoend* B) 

SYS=1ST 

SYSIST. 

(1 in lower half 
of RJ word 

& SYSERR* (Append* 

(other than 1 in 
lower half of RJ word) 

1 SYSERR* (Append* 

TAN 

TAN (x) 

X valid^ not an odd mutiole of pi72 216 175 142 109 

S. TAN, (x) 

TAN. 

TAN. (X) 


x = 0 

|x|< 2 **Y^ x=n(pi/2)»y. pi74<x<pi/4 

617 

155 

n=0 

1065 

156 

n odd 

1061 

147 

n even 

1061 

157 


299 106 

B) 

692 377 239 

91 
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RoutIne 

Entry Points 

Arguments 

S,Times at Entry Points (argument) Times for CYBER 

173 72 73 74 76 


TANH 

TANH (x) 

X valid 

I TANH. (X) 


TANH. 

TANH, (x) 

X valid andt 
|x|<.55 
.55<|x|<17.1 
|x|>17.1 

XOR 

XOR (x(l)«..., x(n)) 

n~? 
n = 3 
n 5s 4 
n = 5 

each additional argument 


XTOO* 

XT00$ (x,y> 

( 0 .»«.} 

I SYSERR. (Apoend. B) 
(O.fX)fX valid x>0 
x<0 

& SYSERR. (Aopend. B) 
(x»y)x<0t X valid 

& SYSERR. (Apoend. B) 
x,y valid«x<0« y*Ioqx>741.67 
I OLOG, (x) 

1 SYSERR. (Apoend, 8) 
y»logx<741.67 

% OLOG. (xl 
& EXP, (y*log x» 


XTOO. 

XTon, lx,y) 
x=0. 

(Xyy) valid* 
1 OLOG, 
1 DEXP. 


x#0 

(x) 

(y*(og x) 


812 

970 

388 


98 

75 

65 

58 



15 3 
210 
126 



213 

164 

96 

100 

276 

213 

117 

118 

340 

262 

139 

144 

404 

309 

160 

156 

64 

49 

21 

19 


445 

341 

197 

147 

476 

454 

389 

343 

158 

199 

2 04 
147 

403 

304 

167 

132 

753 

606 

280 

188 

684 

558 

239 

149 


129 

99 

66 

62 

406 

352 

120 

79 
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Routine 

Entry Points 

Arguments 

STimes at Entry Points iargumentl 


Times for CYBER 
173 72 73 74 


XTOI* 

XTOT* lx,nl 

1 XTOI. ffx.nll 


XTOI. 

XTOI. lx.n) 

X valId n vafid* n>Q M^en x^S 
nsO 

n<0*r«place n by -n and 

X by l/x« addt 

n*l 

n=2 

n=4 


76 


188 
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Routine 

Entry Points 

Arguments 


STimes at Entry Points (argument! 

Times 
173 72 

for 

73 

CYBER 

74 

76 

XTOY* 

XTOY* (x,y» 





( 0 • 9 0 • ) 

8, SYSERR. (Apoand* B) 

283 

179 

186 

155 

( 0 * 9 X}vatid«x>0 

396 

330 

92 

198 

(0•« x), X val id»x<0 

1 SYSERR* (Append. B) 

368 

284 

189 

155 

X.y va(id * x< 0 

S SYSERR. (Append. B) 

309 

243 

157 

114 

x.y valid, x>0 valid, x^O 
i ALOG. (xl 

X EXP. (y*log x! 

XTOY. 

XTOY. (x,Z) 

399 

315 

201 

141 

(0,x) valid, x>0 

80 

62 

58 

53 

(x,yl valid,x#0 
i ALOG. (x) 

1. EXP. (v»log X) 

XTOZ* 

XT0Z$ (x, 2 ) 

(0.,Z» 

174 

150 

45 

53 

z valid,re(zl>0 

401 

341 

135 

178 

z valid, re(z1<0 

8. SYSERR. (Append. Bl 

398 

296 

212 

121 

z va1id,re(z)-0 

8 SYSERR. (Apoend. B) 

355 

294 

180 

121 

x,z valid, x<0 

8> SYSERR. (Append. B) 
x,z valid. 

312 

240 

156 

104 

x>0 re(z)»log x >741.67 

8 ALOG. (xl 

8 SYSERR. (Append. 8) 

632 

469 

251 

130 

(x,z} valid, xXO 

8 ALOG. (xl 

8 COS.SIN (iffl(zl*log xl 

8 EXP, (re(zl*log xl 

XTOZ. 

XTOZ. (x,zl 
(0., z) 

705 

573 

221 

85 

z valid, re(zl>0 

82 

341 

58 

55 

(x,z) val id, x>0 

8 ALOG. (Xl 

8 EXP. (re(zl»fog x) 

8 COS.SIN (im(z)*fog xl 

476 

422 

94 

91 
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Routine 

entry Points 

Arguments 


STimes at Entry Points fargument} 

Times 

for ! 

CYBER 



173 72 

73 

74 

76 

ZTOI* 





ZTOI$ <z,n) 





(Q.,0.} 

379 

303 

189 

140 

S, SYSERR, (Apoend. B» 





(otxlf x>0 

232 

178 

137 

111 

(0*«x)tx<0 

369 

287 

182 

61 

1 SYSERR, (Append, B» 





z^0« z«n valid 

204 

181 

113 

99 

1. ZTOI. (z,n) 





ZTOI. 





ZTOI. (z,n) 





(z*nl valid n = 0 

85 

66 

51 

54 

n = 1 

233 

230 

115 

85 

r = Z 

710 

602 

179 

125 

n = 3 

725 

614 

178 

122 

n = -1 

656 

571 

151 

118 

n = -2 

1036 

899 

215 

156 

n = -3 

1101 

955 

214 

160 

If n<0» replace n by -n* and add 
n odd 

374 

337 

46 

32 

n even 

327 

291 

36 

32 

If n>3, t=time a(1)fb(1)•log(2)n<t< 
a(2)«-b(2l*log(2nn, where 
adl = 

477. 

295. 

142.5 

86,9 

b(ll = 

233. 

222. 

36.8 

55.0 

al2) = 

162. 

127. 

87.9 

74.5 

b(2) = 

390. 

337. 

62.3 

28.1 

190 
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Figure 1 - Relative error ir SORT afgorithis 
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Figure ? - Error of algorithm for relative 
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Figure 3 - Error in the algorithm used to approximate 

sinfxl over C'*pi/4« pi/41 189 
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Figure 7 - Graph of relative error in the algorithm used 
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FIGURE k 






















































FIGURE 6 

Graph of error in the algorithm used in ARCSTN. for arcsin 
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FIGURE 9 

Error in the atgorithm used in EXP 
over f-log 2/16* log 2/161 

























































FIGURE 10 

Error of approximation in the series Ctruncated) 
for slnh over (-22*.221 





























































FIGURE 11 

Error in the polynonial approxifliation to tanh* over C-12,*12) 
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